Logo ROOT  
Reference Guide
Fitter.cxx
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Author: L. Moneta Mon Sep 4 17:00:10 2006
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7 * *
8 * *
9 **********************************************************************/
10
11// Implementation file for class Fitter
12
13
14#include "Fit/Fitter.h"
15#include "Fit/Chi2FCN.h"
18#include "Math/Minimizer.h"
21#include "Fit/BasicFCN.h"
22#include "Fit/BinData.h"
23#include "Fit/UnBinData.h"
24#include "Fit/FcnAdapter.h"
25#include "Fit/FitConfig.h"
26#include "Fit/FitResult.h"
27#include "Math/Error.h"
28
29#include <memory>
30
31#include "Math/IParamFunction.h"
32
34
35// #include "TMatrixDSym.h"
36// for debugging
37//#include "TMatrixD.h"
38// #include <iomanip>
39
40namespace ROOT {
41
42 namespace Fit {
43
44// use a static variable to get default minimizer options for error def
45// to see if user has changed it later on. If it has not been changed we set
46// for the likelihood method an error def of 0.5
47// t.b.d : multiply likelihood by 2 so have same error def definition as chi2
49
50
51Fitter::Fitter(const std::shared_ptr<FitResult> & result) :
52 fResult(result)
53{
54 if (result->fFitFunc) SetFunction(*fResult->fFitFunc); // this will create also the configuration
55 if (result->fObjFunc) fObjFunction = fResult->fObjFunc;
56 if (result->fFitData) fData = fResult->fFitData;
57}
58
59void Fitter::SetFunction(const IModelFunction & func, bool useGradient)
60{
61
62 fUseGradient = useGradient;
63 if (fUseGradient) {
64 const IGradModelFunction * gradFunc = dynamic_cast<const IGradModelFunction*>(&func);
65 if (gradFunc) {
66 SetFunction(*gradFunc, true);
67 return;
68 }
69 else {
70 MATH_WARN_MSG("Fitter::SetFunction","Requested function does not provide gradient - use it as non-gradient function ");
71 }
72 }
73 fUseGradient = false;
74
75 // set the fit model function (clone the given one and keep a copy )
76 //std::cout << "set a non-grad function" << std::endl;
77
78 fFunc = std::shared_ptr<IModelFunction>(dynamic_cast<IModelFunction *>(func.Clone() ) );
79 assert(fFunc);
80
81 // creates the parameter settings
83 fFunc_v.reset();
84}
85
86void Fitter::SetFunction(const IModel1DFunction & func, bool useGradient)
87{
88 fUseGradient = useGradient;
89 if (fUseGradient) {
90 const IGradModel1DFunction * gradFunc = dynamic_cast<const IGradModel1DFunction*>(&func);
91 if (gradFunc) {
92 SetFunction(*gradFunc, true);
93 return;
94 }
95 else {
96 MATH_WARN_MSG("Fitter::SetFunction","Requested function does not provide gradient - use it as non-gradient function ");
97 }
98 }
99 fUseGradient = false;
100 //std::cout << "set a 1d function" << std::endl;
101
102 // function is cloned when creating the adapter
103 fFunc = std::shared_ptr<IModelFunction>(new ROOT::Math::MultiDimParamFunctionAdapter(func));
104
105 // creates the parameter settings
107 fFunc_v.reset();
108}
109
110void Fitter::SetFunction(const IGradModelFunction & func, bool useGradient)
111{
112 fUseGradient = useGradient;
113 //std::cout << "set a grad function" << std::endl;
114 // set the fit model function (clone the given one and keep a copy )
115 fFunc = std::shared_ptr<IModelFunction>( dynamic_cast<IGradModelFunction *> ( func.Clone() ) );
116 assert(fFunc);
117
118 // creates the parameter settings
120 fFunc_v.reset();
121}
122
123
124void Fitter::SetFunction(const IGradModel1DFunction & func, bool useGradient)
125{
126 //std::cout << "set a 1d grad function" << std::endl;
127 fUseGradient = useGradient;
128 // function is cloned when creating the adapter
129 fFunc = std::shared_ptr<IModelFunction>(new ROOT::Math::MultiDimParamGradFunctionAdapter(func));
130
131 // creates the parameter settings
133 fFunc_v.reset();
134}
135
136
137bool Fitter::DoSetFCN(bool extFcn, const ROOT::Math::IMultiGenFunction & fcn, const double * params, unsigned int dataSize, bool chi2fit) {
138 // Set the objective function for the fit. First parameter specifies if function object is managed external or internal.
139 // In case of an internal function object we need to clone because it is a temporary one
140 // if params is not NULL create the parameter settings
141 fUseGradient = false;
142 unsigned int npar = fcn.NDim();
143 if (npar == 0) {
144 MATH_ERROR_MSG("Fitter::SetFCN","FCN function has zero parameters ");
145 return false;
146 }
147 if (params != 0 )
148 fConfig.SetParamsSettings(npar, params);
149 else {
150 if ( fConfig.ParamsSettings().size() != npar) {
151 MATH_ERROR_MSG("Fitter::SetFCN","wrong fit parameter settings");
152 return false;
153 }
154 }
155
156 fBinFit = chi2fit;
157 fDataSize = dataSize;
158
159 // store external provided FCN without cloning it
160 // it will be cloned in fObjFunc after the fit
161 if (extFcn) {
162 fExtObjFunction = &fcn;
163 fObjFunction.reset();
164 }
165 else {
166 // case FCN is built from Minuit interface so function object is created internally in Fitter class
167 // and needs to be cloned and managed
168 fExtObjFunction = nullptr;
169 fObjFunction.reset(fcn.Clone());
170 }
171
172 // in case a model function and data exists from a previous fit - reset shared-ptr
173 if (fResult && fResult->FittedFunction() == 0 && fFunc) fFunc.reset();
174 if (fData) fData.reset();
175
176 return true;
177}
178bool Fitter::SetFCN(const ROOT::Math::IMultiGenFunction & fcn, const double * params, unsigned int dataSize, bool chi2fit) {
179 // set the objective function for the fit
180 return DoSetFCN(true, fcn, params, dataSize, chi2fit);
181}
182bool Fitter::SetFCN(const ROOT::Math::IMultiGenFunction &fcn, const IModelFunction & func, const double *params, unsigned int dataSize, bool chi2fit) {
183 // set the objective function for the fit and a model function
184 if (!SetFCN(fcn, params, dataSize, chi2fit) ) return false;
185 // need to set fFunc afterwards because SetFCN could reset fFUnc
186 fFunc = std::shared_ptr<IModelFunction>(dynamic_cast<IModelFunction *>(func.Clone()));
187 return (fFunc != nullptr);
188}
189
190bool Fitter::SetFCN(const ROOT::Math::IMultiGradFunction &fcn, const double *params, unsigned int dataSize,
191 bool chi2fit)
192{
193 // set the objective function for the fit
194 // if params is not NULL create the parameter settings
195 if (!SetFCN(static_cast<const ROOT::Math::IMultiGenFunction &>(fcn), params, dataSize, chi2fit))
196 return false;
197 fUseGradient = true;
198 return true;
199}
200
201bool Fitter::SetFCN(const ROOT::Math::IMultiGradFunction &fcn, const IModelFunction &func, const double *params,
202 unsigned int dataSize, bool chi2fit)
203{
204 // set the objective function for the fit and a model function
205 if (!SetFCN(fcn, params, dataSize, chi2fit) ) return false;
206 fFunc = std::shared_ptr<IModelFunction>(dynamic_cast<IModelFunction *>(func.Clone()));
207 return (fFunc != nullptr);
208}
209
210bool Fitter::SetFCN(const ROOT::Math::FitMethodFunction &fcn, const double *params)
211{
212 // set the objective function for the fit
213 // if params is not NULL create the parameter settings
214 bool chi2fit = (fcn.Type() == ROOT::Math::FitMethodFunction::kLeastSquare);
215 if (!SetFCN(fcn, params, fcn.NPoints(), chi2fit))
216 return false;
217 fUseGradient = false;
218 fFitType = fcn.Type();
219 return true;
220}
221
222bool Fitter::SetFCN(const ROOT::Math::FitMethodGradFunction &fcn, const double *params)
223{
224 // set the objective function for the fit
225 // if params is not NULL create the parameter settings
227 if (!SetFCN(fcn, params, fcn.NPoints(), chi2fit))
228 return false;
229 fUseGradient = true;
230 fFitType = fcn.Type();
231 return true;
232}
233
234bool Fitter::FitFCN(const BaseFunc &fcn, const double *params, unsigned int dataSize, bool chi2fit)
235{
236 // fit a user provided FCN function
237 // create fit parameter settings
238 if (!SetFCN(fcn, params, dataSize, chi2fit))
239 return false;
240 return FitFCN();
241}
242
243bool Fitter::FitFCN(const BaseGradFunc &fcn, const double *params, unsigned int dataSize, bool chi2fit)
244{
245 // fit a user provided FCN gradient function
246
247 if (!SetFCN(fcn, params, dataSize, chi2fit))
248 return false;
249 return FitFCN();
250}
251
252bool Fitter::FitFCN(const ROOT::Math::FitMethodFunction &fcn, const double *params)
253{
254 // fit using the passed objective function for the fit
255 if (!SetFCN(fcn, params))
256 return false;
257 return FitFCN();
258}
259
260bool Fitter::FitFCN(const ROOT::Math::FitMethodGradFunction &fcn, const double *params)
261{
262 // fit using the passed objective function for the fit
263 if (!SetFCN(fcn, params))
264 return false;
265 return FitFCN();
266}
267
268bool Fitter::SetFCN(MinuitFCN_t fcn, int npar, const double *params, unsigned int dataSize, bool chi2fit)
269{
270 // set TMinuit style FCN type (global function pointer)
271 // create corresponfing objective function from that function
272
273 if (npar == 0) {
274 npar = fConfig.ParamsSettings().size();
275 if (npar == 0) {
276 MATH_ERROR_MSG("Fitter::FitFCN", "Fit Parameter settings have not been created ");
277 return false;
278 }
279 }
280
281 ROOT::Fit::FcnAdapter newFcn(fcn, npar);
282 return DoSetFCN(false,newFcn, params, dataSize, chi2fit);
283}
284
285bool Fitter::FitFCN(MinuitFCN_t fcn, int npar, const double *params, unsigned int dataSize, bool chi2fit)
286{
287 // fit using Minuit style FCN type (global function pointer)
288 // create corresponfing objective function from that function
289 if (!SetFCN(fcn, npar, params, dataSize, chi2fit))
290 return false;
291 fUseGradient = false;
292 return FitFCN();
293}
294
296{
297 // fit using the previously set FCN function
298
299
300 if (!fExtObjFunction && !fObjFunction) {
301 MATH_ERROR_MSG("Fitter::FitFCN", "Objective function has not been set");
302 return false;
303 }
304 // look if FCN is of a known type and we can get retrieve the model function and data objects
305 if (!fFunc || !fData)
306 ExamineFCN();
307 // init the minimizer
308 if (!DoInitMinimizer())
309 return false;
310 // perform the minimization
311 return DoMinimization();
312}
313
315{
316 // evaluate the FCN using the stored values in fConfig
317
318 if (fFunc && fResult->FittedFunction() == 0)
319 fFunc.reset();
320
321 if (!ObjFunction()) {
322 MATH_ERROR_MSG("Fitter::FitFCN", "Objective function has not been set");
323 return false;
324 }
325 // create a Fit result from the fit configuration
326 fResult = std::make_shared<ROOT::Fit::FitResult>(fConfig);
327 // evaluate one time the FCN
328 double fcnval = (*ObjFunction())(fResult->GetParams());
329 // update fit result
330 fResult->fVal = fcnval;
331 fResult->fNCalls++;
332 return true;
333}
334
336{
337
338 // perform a chi2 fit on a set of binned data
339 std::shared_ptr<BinData> data = std::dynamic_pointer_cast<BinData>(fData);
340 assert(data);
341
342 // check function
343 if (!fFunc && !fFunc_v) {
344 MATH_ERROR_MSG("Fitter::DoLeastSquareFit", "model function is not set");
345 return false;
346 } else {
347
348#ifdef DEBUG
349 std::cout << "Fitter ParamSettings " << Config().ParamsSettings()[3].IsBound() << " lower limit "
350 << Config().ParamsSettings()[3].LowerLimit() << " upper limit "
351 << Config().ParamsSettings()[3].UpperLimit() << std::endl;
352#endif
353
354 fBinFit = true;
355 fDataSize = data->Size();
356 // check if fFunc provides gradient
357 if (!fUseGradient) {
358 // do minimization without using the gradient
359 if (fFunc_v) {
360 return DoMinimization(std::make_unique<Chi2FCN<BaseFunc, IModelFunction_v>>(data, fFunc_v, executionPolicy));
361 } else {
362 return DoMinimization(std::make_unique<Chi2FCN<BaseFunc>>(data, fFunc, executionPolicy));
363 }
364 } else {
365 // use gradient
367 MATH_INFO_MSG("Fitter::DoLeastSquareFit", "use gradient from model function");
368
369 if (fFunc_v) {
370 std::shared_ptr<IGradModelFunction_v> gradFun = std::dynamic_pointer_cast<IGradModelFunction_v>(fFunc_v);
371 if (gradFun) {
372 return DoMinimization(std::make_unique<Chi2FCN<BaseGradFunc, IModelFunction_v>>(data, gradFun, executionPolicy));
373 }
374 } else {
375 std::shared_ptr<IGradModelFunction> gradFun = std::dynamic_pointer_cast<IGradModelFunction>(fFunc);
376 if (gradFun) {
377 return DoMinimization(std::make_unique<Chi2FCN<BaseGradFunc>>(data, gradFun, executionPolicy));
378 }
379 }
380 MATH_ERROR_MSG("Fitter::DoLeastSquareFit", "wrong type of function - it does not provide gradient");
381 }
382 }
383 return false;
384}
385
386bool Fitter::DoBinnedLikelihoodFit(bool extended, const ROOT::EExecutionPolicy &executionPolicy)
387{
388 // perform a likelihood fit on a set of binned data
389 // The fit is extended (Poisson logl_ by default
390
391 std::shared_ptr<BinData> data = std::dynamic_pointer_cast<BinData>(fData);
392 assert(data);
393
394 bool useWeight = fConfig.UseWeightCorrection();
395
396 // check function
397 if (!fFunc && !fFunc_v) {
398 MATH_ERROR_MSG("Fitter::DoBinnedLikelihoodFit", "model function is not set");
399 return false;
400 }
401
402 // logl fit (error should be 0.5) set if different than default values (of 1)
405 }
406
407 if (useWeight && fConfig.MinosErrors()) {
408 MATH_INFO_MSG("Fitter::DoBinnedLikelihoodFit", "MINOS errors cannot be computed in weighted likelihood fits");
409 fConfig.SetMinosErrors(false);
410 }
411
412 fBinFit = true;
413 fDataSize = data->Size();
414
415 if (!fUseGradient) {
416 // do minimization without using the gradient
417 if (fFunc_v) {
418 // create a chi2 function to be used for the equivalent chi-square
420 auto logl = std::make_unique<PoissonLikelihoodFCN<BaseFunc, IModelFunction_v>>(data, fFunc_v, useWeight, extended, executionPolicy);
421 return (useWeight) ? DoWeightMinimization(std::move(logl),&chi2) : DoMinimization(std::move(logl),&chi2);
422 } else {
423 // create a chi2 function to be used for the equivalent chi-square
425 auto logl = std::make_unique<PoissonLikelihoodFCN<BaseFunc>>(data, fFunc, useWeight, extended, executionPolicy);
426 return (useWeight) ? DoWeightMinimization(std::move(logl),&chi2) : DoMinimization(std::move(logl),&chi2);
427 }
428 } else {
430 MATH_INFO_MSG("Fitter::DoLikelihoodFit", "use gradient from model function");
431 // not-extended is not implemented in this case
432 if (!extended) {
433 MATH_WARN_MSG("Fitter::DoBinnedLikelihoodFit",
434 "Not-extended binned fit with gradient not yet supported - do an extended fit");
435 extended = true;
436 }
437 if (fFunc_v) {
438 // create a chi2 function to be used for the equivalent chi-square
440 std::shared_ptr<IGradModelFunction_v> gradFun = std::dynamic_pointer_cast<IGradModelFunction_v>(fFunc_v);
441 if (!gradFun) {
442 MATH_ERROR_MSG("Fitter::DoBinnedLikelihoodFit", "wrong type of function - it does not provide gradient");
443 return false;
444 }
445 auto logl = std::make_unique<PoissonLikelihoodFCN<BaseGradFunc, IModelFunction_v>>(data, gradFun, useWeight, extended, executionPolicy);
446 // do minimization
447 return (useWeight) ? DoWeightMinimization(std::move(logl),&chi2) : DoMinimization(std::move(logl),&chi2);
448 } else {
449 // create a chi2 function to be used for the equivalent chi-square
451 // check if fFunc provides gradient
452 std::shared_ptr<IGradModelFunction> gradFun = std::dynamic_pointer_cast<IGradModelFunction>(fFunc);
453 if (!gradFun) {
454 MATH_ERROR_MSG("Fitter::DoBinnedLikelihoodFit", "wrong type of function - it does not provide gradient");
455 return false;
456 }
457 // use gradient for minimization
458 auto logl = std::make_unique<PoissonLikelihoodFCN<BaseGradFunc>>(data, gradFun, useWeight, extended, executionPolicy);
459 // do minimization
460 return (useWeight) ? DoWeightMinimization(std::move(logl),&chi2) : DoMinimization(std::move(logl),&chi2);
461 }
462 }
463 return false;
464}
465
466bool Fitter::DoUnbinnedLikelihoodFit(bool extended, const ROOT::EExecutionPolicy &executionPolicy) {
467 // perform a likelihood fit on a set of unbinned data
468
469 std::shared_ptr<UnBinData> data = std::dynamic_pointer_cast<UnBinData>(fData);
470 assert(data);
471
472 bool useWeight = fConfig.UseWeightCorrection();
473
474 if (!fFunc && !fFunc_v) {
475 MATH_ERROR_MSG("Fitter::DoUnbinnedLikelihoodFit","model function is not set");
476 return false;
477 }
478
479 if (useWeight && fConfig.MinosErrors() ) {
480 MATH_INFO_MSG("Fitter::DoUnbinnedLikelihoodFit","MINOS errors cannot be computed in weighted likelihood fits");
481 fConfig.SetMinosErrors(false);
482 }
483
484
485 fBinFit = false;
486 fDataSize = data->Size();
487
488#ifdef DEBUG
489 int ipar = 0;
490 std::cout << "Fitter ParamSettings " << Config().ParamsSettings()[ipar].IsBound() << " lower limit " << Config().ParamsSettings()[ipar].LowerLimit() << " upper limit " << Config().ParamsSettings()[ipar].UpperLimit() << std::endl;
491#endif
492
493 // logl fit (error should be 0.5) set if different than default values (of 1)
496 }
497
498 if (!fUseGradient) {
499 // do minimization without using the gradient
500 if (fFunc_v ){
501 auto logl = std::make_unique<LogLikelihoodFCN<BaseFunc, IModelFunction_v>>(data, fFunc_v, useWeight, extended, executionPolicy);
502 // do minimization
503 return (useWeight) ? DoWeightMinimization(std::move(logl)) : DoMinimization(std::move(logl));
504 } else {
505 auto logl = std::make_unique<LogLikelihoodFCN<BaseFunc>>(data, fFunc, useWeight, extended, executionPolicy);
506 return (useWeight) ? DoWeightMinimization(std::move(logl)) : DoMinimization(std::move(logl));
507 }
508 } else {
509 // use gradient : check if fFunc provides gradient
511 MATH_INFO_MSG("Fitter::DoUnbinnedLikelihoodFit", "use gradient from model function");
512 if (extended) {
513 MATH_WARN_MSG("Fitter::DoUnbinnedLikelihoodFit",
514 "Extended unbinned fit with gradient not yet supported - do a not-extended fit");
515 extended = false;
516 }
517 if (fFunc_v) {
518 std::shared_ptr<IGradModelFunction_v> gradFun = std::dynamic_pointer_cast<IGradModelFunction_v>(fFunc_v);
519 if (!gradFun) {
520 MATH_ERROR_MSG("Fitter::DoUnbinnedLikelihoodFit", "wrong type of function - it does not provide gradient");
521 return false;
522 }
523 auto logl = std::make_unique<LogLikelihoodFCN<BaseGradFunc, IModelFunction_v>>(data, gradFun, useWeight, extended, executionPolicy);
524 return (useWeight) ? DoWeightMinimization(std::move(logl)) : DoMinimization(std::move(logl));
525 } else {
526 std::shared_ptr<IGradModelFunction> gradFun = std::dynamic_pointer_cast<IGradModelFunction>(fFunc);
527 if (!gradFun) {
528 MATH_ERROR_MSG("Fitter::DoUnbinnedLikelihoodFit", "wrong type of function - it does not provide gradient");
529 return false;
530 }
531 auto logl = std::make_unique<LogLikelihoodFCN<BaseGradFunc>>(data, gradFun, useWeight, extended, executionPolicy);
532 return (useWeight) ? DoWeightMinimization(std::move(logl)) : DoMinimization(std::move(logl));
533 }
534 }
535 return false;
536}
537
538
540
541 std::shared_ptr<BinData> data = std::dynamic_pointer_cast<BinData>(fData);
542 assert(data);
543
544 // perform a linear fit on a set of binned data
545 std::string prevminimizer = fConfig.MinimizerType();
546 fConfig.SetMinimizer("Linear");
547
548 fBinFit = true;
549
550 bool ret = DoLeastSquareFit();
551 fConfig.SetMinimizer(prevminimizer.c_str());
552 return ret;
553}
554
555
557 // compute the Hesse errors according to configuration
558 // set in the parameters and append value in fit result
559 if (!ObjFunction()) {
560 MATH_ERROR_MSG("Fitter::CalculateHessErrors","Objective function has not been set");
561 return false;
562 }
563
564 // need a special treatment in case of weighted likelihood fit
565 // (not yet implemented)
566 if (fFitType == 2 && fConfig.UseWeightCorrection() ) {
567 MATH_ERROR_MSG("Fitter::CalculateHessErrors","Re-computation of Hesse errors not implemented for weighted likelihood fits");
568 MATH_INFO_MSG("Fitter::CalculateHessErrors","Do the Fit using configure option FitConfig::SetParabErrors()");
569 return false;
570 }
571
572 // a fit Result pointer must exist when a minimizer exists
573 if (fMinimizer && !fResult ) {
574 MATH_ERROR_MSG("Fitter::CalculateHessErrors", "FitResult has not been created");
575 return false;
576 }
577
578 // update minimizer (recreate if not done or if name has changed
580 MATH_ERROR_MSG("Fitter::CalculateHessErrors","Error re-initializing the minimizer");
581 return false;
582 }
583
584 if (!fMinimizer ) {
585 // this should not happen
586 MATH_ERROR_MSG("Fitter::CalculateHessErrors", "Need to do a fit before calculating the errors");
587 assert(false);
588 return false;
589 }
590
591 //run Hesse
592 bool ret = fMinimizer->Hesse();
593 if (!ret) MATH_WARN_MSG("Fitter::CalculateHessErrors","Error when calculating Hessian");
594
595 // update minimizer results with what comes out from Hesse
596 // in case is empty - create from a FitConfig
597 if (fResult->IsEmpty() )
599
600 // update obj function in case it was an external one
602 fResult->fObjFunc = fObjFunction;
603
604 // re-give a minimizer instance in case it has been changed
605 ret |= fResult->Update(fMinimizer, fConfig, ret);
606
607 // when possible get ncalls from FCN and set in fit result
609 fResult->fNCalls = GetNCallsFromFCN();
610 }
611
612 // set also new errors in FitConfig
614
615 return ret;
616}
617
618
620 // compute the Minos errors according to configuration
621 // set in the parameters and append value in fit result
622 // normally Minos errors are computed just after the minimization
623 // (in DoMinimization) aftewr minimizing if the
624 // FitConfig::MinosErrors() flag is set
625
626 if (!fMinimizer) {
627 MATH_ERROR_MSG("Fitter::CalculateMinosErrors","Minimizer does not exist - cannot calculate Minos errors");
628 return false;
629 }
630
631 if (!fResult || fResult->IsEmpty() ) {
632 MATH_ERROR_MSG("Fitter::CalculateMinosErrors","Invalid Fit Result - cannot calculate Minos errors");
633 return false;
634 }
635
636 if (fFitType == 2 && fConfig.UseWeightCorrection() ) {
637 MATH_ERROR_MSG("Fitter::CalculateMinosErrors","Computation of MINOS errors not implemented for weighted likelihood fits");
638 return false;
639 }
640
641 // update minimizer (but cannot re-create in this case). Must use an existing one
642 if (!DoUpdateMinimizerOptions(false)) {
643 MATH_ERROR_MSG("Fitter::CalculateHessErrors","Error re-initializing the minimizer");
644 return false;
645 }
646
647 // set flag to compute Minos error to false in FitConfig to avoid that
648 // following minimizaiton calls perform unwanted Minos error calculations
649 /// fConfig.SetMinosErrors(false);
650
651
652 const std::vector<unsigned int> & ipars = fConfig.MinosParams();
653 unsigned int n = (ipars.size() > 0) ? ipars.size() : fResult->Parameters().size();
654 bool ok = false;
655
656 int iparNewMin = 0;
657 int iparMax = n;
658 int iter = 0;
659 // rerun minos for the parameters run before a new Minimum has been found
660 do {
661 if (iparNewMin > 0)
662 MATH_INFO_MSG("Fitter::CalculateMinosErrors","Run again Minos for some parameters because a new Minimum has been found");
663 iparNewMin = 0;
664 for (int i = 0; i < iparMax; ++i) {
665 double elow, eup;
666 unsigned int index = (ipars.size() > 0) ? ipars[i] : i;
667 bool ret = fMinimizer->GetMinosError(index, elow, eup);
668 // flags case when a new minimum has been found
669 if ((fMinimizer->MinosStatus() & 8) != 0) {
670 iparNewMin = i;
671 }
672 if (ret)
673 fResult->SetMinosError(index, elow, eup);
674 ok |= ret;
675 }
676
677 iparMax = iparNewMin;
678 iter++; // to avoid infinite looping
679 }
680 while( iparNewMin > 0 && iter < 10);
681 if (!ok) {
682 MATH_ERROR_MSG("Fitter::CalculateMinosErrors","Minos error calculation failed for all the selected parameters");
683 }
684
685 // update obj function in case it was an external one
687 fResult->fObjFunc = fObjFunction;
688
689 // re-give a minimizer instance in case it has been changed
690 // but maintain previous valid status. Do not set result to false if minos failed
691 ok &= fResult->Update(fMinimizer, fConfig, fResult->IsValid());
692
693 return ok;
694}
695
696
697
698// traits for distinguishing fit methods functions from generic objective functions
699template<class Func>
701 static unsigned int NCalls(const Func & ) { return 0; }
702 static int Type(const Func & ) { return -1; }
703 static bool IsGrad() { return false; }
704};
705template<>
707 static unsigned int NCalls(const ROOT::Math::FitMethodFunction & f ) { return f.NCalls(); }
708 static int Type(const ROOT::Math::FitMethodFunction & f) { return f.Type(); }
709 static bool IsGrad() { return false; }
710};
711template<>
713 static unsigned int NCalls(const ROOT::Math::FitMethodGradFunction & f ) { return f.NCalls(); }
714 static int Type(const ROOT::Math::FitMethodGradFunction & f) { return f.Type(); }
715 static bool IsGrad() { return true; }
716};
717
719 //initialize minimizer by creating it
720 // and set there the objective function
721 // obj function must have been set before
722 auto objFunction = ObjFunction();
723 if (!objFunction) {
724 MATH_ERROR_MSG("Fitter::DoInitMinimizer","Objective function has not been set");
725 return false;
726 }
727
728 // check configuration and objective function
729 if ( fConfig.ParamsSettings().size() != objFunction->NDim() ) {
730 MATH_ERROR_MSG("Fitter::DoInitMinimizer","wrong function dimension or wrong size for FitConfig");
731 return false;
732 }
733
734 // create first Minimizer
735 // using an auto_Ptr will delete the previous existing one
736 fMinimizer = std::shared_ptr<ROOT::Math::Minimizer> ( fConfig.CreateMinimizer() );
737 if (fMinimizer.get() == 0) {
738 MATH_ERROR_MSG("Fitter::DoInitMinimizer","Minimizer cannot be created");
739 return false;
740 }
741
742 // in case of gradient function one needs to downcast the pointer
743 if (fUseGradient) {
744 const ROOT::Math::IMultiGradFunction * gradfcn = dynamic_cast<const ROOT::Math::IMultiGradFunction *> (objFunction );
745 if (!gradfcn) {
746 MATH_ERROR_MSG("Fitter::DoInitMinimizer","wrong type of function - it does not provide gradient");
747 return false;
748 }
749 fMinimizer->SetFunction( *gradfcn);
750 }
751 else
752 fMinimizer->SetFunction( *objFunction);
753
754
755 fMinimizer->SetVariables(fConfig.ParamsSettings().begin(), fConfig.ParamsSettings().end() );
756
757 // if requested parabolic error do correct error analysis by the minimizer (call HESSE)
758 if (fConfig.ParabErrors()) fMinimizer->SetValidError(true);
759
760 return true;
761
762}
763
764bool Fitter::DoUpdateMinimizerOptions(bool canDifferentMinim ) {
765 // update minimizer options when re-doing a Fit or computing Hesse or Minos errors
766
767
768 // create a new minimizer if it is different type
769 // minimizer type string stored in FitResult is "minimizer name" + " / " + minimizer algo
770 std::string newMinimType = fConfig.MinimizerName();
771 if (fMinimizer && fResult && newMinimType != fResult->MinimizerType()) {
772 // if a different minimizer is allowed (e.g. when calling Hesse)
773 if (canDifferentMinim) {
774 std::string msg = "Using now " + newMinimType;
775 MATH_INFO_MSG("Fitter::DoUpdateMinimizerOptions: ", msg.c_str());
776 if (!DoInitMinimizer() )
777 return false;
778 }
779 else {
780 std::string msg = "Cannot change minimizer. Continue using " + fResult->MinimizerType();
781 MATH_WARN_MSG("Fitter::DoUpdateMinimizerOptions",msg.c_str());
782 }
783 }
784
785 // create minimizer if it was not done before
786 if (!fMinimizer) {
787 if (!DoInitMinimizer())
788 return false;
789 }
790
791 // set new minimizer options (but not functions and parameters)
792 fMinimizer->SetOptions(fConfig.MinimizerOptions());
793 return true;
794}
795
797 // perform the minimization (assume we have already initialized the minimizer)
798
799 assert(fMinimizer );
800
801 bool isValid = fMinimizer->Minimize();
802
803 if (!fResult) fResult = std::make_shared<FitResult>();
804
805 fResult->FillResult(fMinimizer,fConfig, fFunc, isValid, fDataSize, fBinFit, chi2func );
806
807 // if requested run Minos after minimization
808 if (isValid && fConfig.MinosErrors()) {
809 // minos error calculation will update also FitResult
811 }
812
813 // when possible get number of calls from FCN and set in fit result
815 fResult->fNCalls = GetNCallsFromFCN();
816 }
817
818 // fill information in fit result
819 // if using an external obj function clone it for storing in FitResult
821 fResult->fObjFunc = fObjFunction;
822 fResult->fFitData = fData;
823
824#ifdef DEBUG
825 std::cout << "ROOT::Fit::Fitter::DoMinimization : ncalls = " << fResult->fNCalls << " type of objfunc " << fFitFitResType << " typeid: " << typeid(*fObjFunction).name() << " use gradient " << fUseGradient << std::endl;
826#endif
827
829 fResult->NormalizeErrors();
830
831 // set also new parameter values and errors in FitConfig
832 if (fConfig.UpdateAfterFit() && isValid) DoUpdateFitConfig();
833
834 return isValid;
835}
836template<class ObjFunc_t>
837bool Fitter::DoMinimization(std::unique_ptr<ObjFunc_t> objFunc, const ROOT::Math::IMultiGenFunction * chi2func) {
838 // perform the minimization initializing the minimizer starting from a given obj function
839 fFitType = objFunc->Type();
840 fExtObjFunction = nullptr;
841 fObjFunction = std::move(objFunc);
842 if (!DoInitMinimizer()) return false;
843 return DoMinimization(chi2func);
844}
845template<class ObjFunc_t>
846bool Fitter::DoWeightMinimization(std::unique_ptr<ObjFunc_t> objFunc, const ROOT::Math::IMultiGenFunction * chi2func) {
847 // perform the minimization initializing the minimizer starting from a given obj function
848 // and apply afterwards the correction for weights. This applyies only for logL fitting
849 this->fFitType = objFunc->Type();
850 fExtObjFunction = nullptr;
851 fObjFunction = std::move(objFunc);
852 if (!DoInitMinimizer()) return false;
853 if (!DoMinimization(chi2func)) return false;
854 objFunc->UseSumOfWeightSquare();
855 return ApplyWeightCorrection(*objFunc);
856}
857
858
860 // update the fit configuration after a fit using the obtained result
861 if (fResult->IsEmpty() || !fResult->IsValid() ) return;
862 for (unsigned int i = 0; i < fConfig.NPar(); ++i) {
864 par.SetValue( fResult->Value(i) );
865 if (fResult->Error(i) > 0) par.SetStepSize( fResult->Error(i) );
866 }
867}
868
870 // retrieve ncalls from the fit method functions
871 // this function is called when minimizer does not provide a way of returning the nnumber of function calls
872 int ncalls = 0;
873 if (!fUseGradient) {
874 const ROOT::Math::FitMethodFunction * fcn = dynamic_cast<const ROOT::Math::FitMethodFunction *>(fObjFunction.get());
875 if (fcn) ncalls = fcn->NCalls();
876 }
877 else {
879 if (fcn) ncalls = fcn->NCalls();
880 }
881 return ncalls;
882}
883
884
885bool Fitter::ApplyWeightCorrection(const ROOT::Math::IMultiGenFunction & loglw2, bool minimizeW2L) {
886 // apply correction for weight square
887 // Compute Hessian of the loglikelihood function using the sum of the weight squared
888 // This method assumes:
889 // - a fit has been done before and a covariance matrix exists
890 // - the objective function is a likelihood function and Likelihood::UseSumOfWeightSquare()
891 // has been called before
892
893 if (fMinimizer.get() == 0) {
894 MATH_ERROR_MSG("Fitter::ApplyWeightCorrection","Must perform first a fit before applying the correction");
895 return false;
896 }
897
898 unsigned int n = loglw2.NDim();
899 // correct errors for weight squared
900 std::vector<double> cov(n*n);
901 bool ret = fMinimizer->GetCovMatrix(&cov[0] );
902 if (!ret) {
903 MATH_ERROR_MSG("Fitter::ApplyWeightCorrection","Previous fit has no valid Covariance matrix");
904 return false;
905 }
906 // need to use new obj function computed with weight-square
907 std::shared_ptr<ROOT::Math::IMultiGenFunction> objFunc(loglw2.Clone());
908 fObjFunction.swap( objFunc );
909
910 // need to re-initialize the minimizer for the changes applied in the
911 // objective functions
912 if (!DoInitMinimizer()) return false;
913
914 //std::cout << "Running Hesse ..." << std::endl;
915
916 // run eventually before a minimization
917 // ignore its error
918 if (minimizeW2L) fMinimizer->Minimize();
919 // run Hesse on the log-likelihood build using sum of weight squared
920 ret = fMinimizer->Hesse();
921 if (!ret) {
922 MATH_ERROR_MSG("Fitter::ApplyWeightCorrection","Error running Hesse on weight2 likelihood - cannot compute errors");
923 return false;
924 }
925
926 if (fMinimizer->CovMatrixStatus() != 3) {
927 MATH_WARN_MSG("Fitter::ApplyWeightCorrection","Covariance matrix for weighted likelihood is not accurate, the errors may be not reliable");
928 if (fMinimizer->CovMatrixStatus() == 2)
929 MATH_WARN_MSG("Fitter::ApplyWeightCorrection","Covariance matrix for weighted likelihood was forced to be defined positive");
930 if (fMinimizer->CovMatrixStatus() <= 0)
931 // probably should have failed before
932 MATH_ERROR_MSG("Fitter::ApplyWeightCorrection","Covariance matrix for weighted likelihood is not valid !");
933 }
934
935 // get Hessian matrix from weight-square likelihood
936 std::vector<double> hes(n*n);
937 ret = fMinimizer->GetHessianMatrix(&hes[0] );
938 if (!ret) {
939 MATH_ERROR_MSG("Fitter::ApplyWeightCorrection","Error retrieving Hesse on weight2 likelihood - cannot compute errors");
940 return false;
941 }
942
943
944 // perform product of matrix cov * hes * cov
945 // since we do not want to add matrix dependence do product by hand
946 // first do hes * cov
947 std::vector<double> tmp(n*n);
948 for (unsigned int i = 0; i < n; ++i) {
949 for (unsigned int j = 0; j < n; ++j) {
950 for (unsigned int k = 0; k < n; ++k)
951 tmp[i*n+j] += hes[i*n + k] * cov[k*n + j];
952 }
953 }
954 // do multiplication now cov * tmp save result
955 std::vector<double> newCov(n*n);
956 for (unsigned int i = 0; i < n; ++i) {
957 for (unsigned int j = 0; j < n; ++j) {
958 for (unsigned int k = 0; k < n; ++k)
959 newCov[i*n+j] += cov[i*n + k] * tmp[k*n + j];
960 }
961 }
962 // update fit result with new corrected covariance matrix
963 unsigned int k = 0;
964 for (unsigned int i = 0; i < n; ++i) {
965 fResult->fErrors[i] = std::sqrt( newCov[i*(n+1)] );
966 for (unsigned int j = 0; j <= i; ++j)
967 fResult->fCovMatrix[k++] = newCov[i *n + j];
968 }
969
970 // restore previous used objective function
971 fObjFunction.swap( objFunc );
972
973 return true;
974}
975
976
977
979 // return a pointer to the binned data used in the fit
980 // works only for chi2 or binned likelihood fits
981 // thus when the objective function stored is a Chi2Func or a PoissonLikelihood
982 // This also set the model function correctly if it has not been set
983
986
989
990 //MATH_INFO_MSG("Fitter::ExamineFCN","Objective function is not of a known type - FitData and ModelFunction objects are not available");
991 return;
992}
993
994 } // end namespace Fit
995
996} // end namespace ROOT
#define MATH_INFO_MSG(loc, str)
Pre-processor macro to report messages which can be configured to use ROOT error or simply an std::io...
Definition: Error.h:77
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:83
#define MATH_WARN_MSG(loc, str)
Definition: Error.h:80
#define f(i)
Definition: RSha256.hxx:104
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
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 result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
BasicFCN class: base class for the objective functions used in the fits It has a reference to the dat...
Definition: BasicFCN.h:40
Chi2FCN class for binned fits using the least square methods.
Definition: Chi2FCN.h:46
const std::vector< unsigned int > & MinosParams() const
return vector of parameter indices for which the Minos Error will be computed
Definition: FitConfig.h:220
bool UpdateAfterFit() const
Update configuration after a fit using the FitResult.
Definition: FitConfig.h:213
void SetMinimizer(const char *type, const char *algo=nullptr)
set minimizer type
Definition: FitConfig.h:181
void SetMinosErrors(bool on=true)
set Minos erros computation to be performed after fitting
Definition: FitConfig.h:231
bool NormalizeErrors() const
flag to check if resulting errors are be normalized according to chi2/ndf
Definition: FitConfig.h:204
bool ParabErrors() const
do analysis for parabolic errors
Definition: FitConfig.h:207
unsigned int NPar() const
number of parameters settings
Definition: FitConfig.h:96
void SetParamsSettings(unsigned int npar, const double *params, const double *vstep=nullptr)
set the parameter settings from number of parameters and a vector of values and optionally step value...
Definition: FitConfig.cxx:135
std::string MinimizerName() const
return Minimizer full name (type / algorithm)
Definition: FitConfig.cxx:239
bool UseWeightCorrection() const
Apply Weight correction for error matrix computation.
Definition: FitConfig.h:216
const std::vector< ROOT::Fit::ParameterSettings > & ParamsSettings() const
get the vector of parameter settings (const method)
Definition: FitConfig.h:86
ROOT::Math::Minimizer * CreateMinimizer()
create a new minimizer according to chosen configuration
Definition: FitConfig.cxx:173
void CreateParamsSettings(const ROOT::Math::IParamMultiFunctionTempl< T > &func)
set the parameter settings from a model function.
Definition: FitConfig.h:109
const std::string & MinimizerType() const
return type of minimizer package
Definition: FitConfig.h:189
const ParameterSettings & ParSettings(unsigned int i) const
get the parameter settings for the i-th parameter (const method)
Definition: FitConfig.h:76
ROOT::Math::MinimizerOptions & MinimizerOptions()
access to the minimizer control parameter (non const method)
Definition: FitConfig.h:167
bool MinosErrors() const
do minos errors analysis on the parameters
Definition: FitConfig.h:210
class containing the result of the fit and all the related information (fitted parameter values,...
Definition: FitResult.h:47
bool DoSetFCN(bool useExtFCN, const ROOT::Math::IMultiGenFunction &fcn, const double *params, unsigned int dataSize, bool chi2fit)
Set Objective function.
Definition: Fitter.cxx:137
bool EvalFCN()
Perform a simple FCN evaluation.
Definition: Fitter.cxx:314
const ROOT::Math::IMultiGenFunction * fExtObjFunction
! pointer to an external FCN
Definition: Fitter.h:579
bool FitFCN()
Perform a fit with the previously set FCN function.
Definition: Fitter.cxx:295
void DoUpdateFitConfig()
Definition: Fitter.cxx:859
const ROOT::Math::IMultiGenFunction * ObjFunction() const
return pointer to used objective function if using an external one use external, otherwise use the on...
Definition: Fitter.h:549
bool DoMinimization(std::unique_ptr< ObjFunc_t > f, const ROOT::Math::IMultiGenFunction *chifunc=nullptr)
do minimization
Definition: Fitter.cxx:837
int fDataSize
size of data sets (need for Fumili or LM fitters)
Definition: Fitter.h:563
bool DoUnbinnedLikelihoodFit(bool extended=false, const ROOT::EExecutionPolicy &executionPolicy=ROOT::EExecutionPolicy::kSequential)
un-binned likelihood fit
Definition: Fitter.cxx:466
std::shared_ptr< ROOT::Math::Minimizer > fMinimizer
! pointer to used minimizer
Definition: Fitter.h:573
bool DoWeightMinimization(std::unique_ptr< ObjFunc_t > f, const ROOT::Math::IMultiGenFunction *chifunc=nullptr)
Definition: Fitter.cxx:846
bool DoBinnedLikelihoodFit(bool extended=true, const ROOT::EExecutionPolicy &executionPolicy=ROOT::EExecutionPolicy::kSequential)
binned likelihood fit
Definition: Fitter.cxx:386
int fFitType
type of fit (0 undefined, 1 least square, 2 likelihood)
Definition: Fitter.h:561
std::shared_ptr< ROOT::Fit::FitData > fData
! pointer to the fit data (binned or unbinned data)
Definition: Fitter.h:575
bool fUseGradient
flag to indicate if using gradient or not
Definition: Fitter.h:555
bool fBinFit
flag to indicate if fit is binned in case of false the fit is unbinned or undefined) flag it is used ...
Definition: Fitter.h:557
bool SetFCN(unsigned int npar, Function &fcn, const double *params=0, unsigned int dataSize=0, bool chi2fit=false)
Set a generic FCN function as a C++ callable object implementing double () (const double *) Note that...
Definition: Fitter.h:660
std::shared_ptr< ROOT::Math::IMultiGenFunction > fObjFunction
! pointer to used objective function
Definition: Fitter.h:577
bool ApplyWeightCorrection(const ROOT::Math::IMultiGenFunction &loglw2, bool minimizeW2L=false)
apply correction in the error matrix for the weights for likelihood fits This method can be called on...
Definition: Fitter.cxx:885
void ExamineFCN()
look at the user provided FCN and get data and model function is they derive from ROOT::Fit FCN class...
Definition: Fitter.cxx:978
const FitConfig & Config() const
access to the fit configuration (const method)
Definition: Fitter.h:442
bool DoLeastSquareFit(const ROOT::EExecutionPolicy &executionPolicy=ROOT::EExecutionPolicy::kSequential)
least square fit
Definition: Fitter.cxx:335
std::shared_ptr< IModelFunction_v > fFunc_v
! copy of the fitted function containing on output the fit result
Definition: Fitter.h:567
std::shared_ptr< ROOT::Fit::FitResult > fResult
! pointer to the object containing the result of the fit
Definition: Fitter.h:571
bool GetDataFromFCN()
internal functions to get data set and model function from FCN useful for fits done with customized F...
Definition: Fitter.h:587
bool CalculateMinosErrors()
perform an error analysis on the result using MINOS To be called only after fitting and when a minimi...
Definition: Fitter.cxx:619
bool DoUpdateMinimizerOptions(bool canDifferentMinim=true)
Definition: Fitter.cxx:764
void SetFunction(const IModelFunction &func, bool useGradient=false)
Set the fitted function (model function) from a parametric function interface.
Definition: Fitter.cxx:59
bool CalculateHessErrors()
perform an error analysis on the result using the Hessian Errors are obtained from the inverse of the...
Definition: Fitter.cxx:556
FitConfig fConfig
fitter configuration (options and parameter settings)
Definition: Fitter.h:565
Fitter()
Default constructor.
Definition: Fitter.h:102
std::shared_ptr< IModelFunction > fFunc
! copy of the fitted function containing on output the fit result
Definition: Fitter.h:569
bool DoLinearFit()
linear least square fit
Definition: Fitter.cxx:539
bool DoInitMinimizer()
Definition: Fitter.cxx:718
int GetNCallsFromFCN()
Definition: Fitter.cxx:869
Class, describing value, limits and step size of the parameters Provides functionality also to set/re...
void SetValue(double val)
set the value
void SetStepSize(double err)
set the step size
FitMethodFunction class Interface for objective functions (like chi2 and likelihood used in the fit) ...
virtual Type_t Type() const
return the type of method, override if needed
virtual unsigned int NPoints() const
return the number of data points used in evaluating the function
virtual unsigned int NCalls() const
return the total number of function calls (override if needed)
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
virtual IBaseFunctionMultiDimTempl< T > * Clone() const =0
Clone a function.
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
Interface (abstract class) for multi-dimensional functions providing a gradient calculation.
Definition: IFunction.h:343
Specialized IParamFunction interface (abstract class) for one-dimensional parametric functions It is ...
Interface (abstract class) for parametric gradient multi-dimensional functions providing in addition ...
Interface (abstract class) for parametric one-dimensional gradient functions providing in addition to...
double ErrorDef() const
error definition
int PrintLevel() const
non-static methods for retrieving options
void SetErrorDef(double err)
set error def
MultiDimParamFunctionAdapter class to wrap a one-dimensional parametric function in a multi dimension...
MultiDimParamGradFunctionAdapter class to wrap a one-dimensional parametric gradient function in a mu...
const Int_t n
Definition: legend1.C:16
TFitResultPtr Fit(FitObject *h1, TF1 *f1, Foption_t &option, const ROOT::Math::MinimizerOptions &moption, const char *goption, ROOT::Fit::DataRange &range)
Definition: HFitImpl.cxx:133
double gDefaultErrorDef
Definition: Fitter.cxx:48
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
static int Type(const ROOT::Math::FitMethodFunction &f)
Definition: Fitter.cxx:708
static unsigned int NCalls(const ROOT::Math::FitMethodFunction &f)
Definition: Fitter.cxx:707
static int Type(const ROOT::Math::FitMethodGradFunction &f)
Definition: Fitter.cxx:714
static unsigned int NCalls(const ROOT::Math::FitMethodGradFunction &f)
Definition: Fitter.cxx:713
static bool IsGrad()
Definition: Fitter.cxx:703
static unsigned int NCalls(const Func &)
Definition: Fitter.cxx:701
static int Type(const Func &)
Definition: Fitter.cxx:702