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
825#ifdef DEBUG
826 std::cout << "ROOT::Fit::Fitter::DoMinimization : ncalls = " << fResult->fNCalls << " type of objfunc " << fFitFitResType << " typeid: " << typeid(*fObjFunction).name() << " use gradient " << fUseGradient << std::endl;
827#endif
828
830 fResult->NormalizeErrors();
831
832 // set also new parameter values and errors in FitConfig
833 if (fConfig.UpdateAfterFit() && isValid) DoUpdateFitConfig();
834
835 return isValid;
836}
837template<class ObjFunc_t>
838bool Fitter::DoMinimization(std::unique_ptr<ObjFunc_t> objFunc, const ROOT::Math::IMultiGenFunction * chi2func) {
839 // perform the minimization initializing the minimizer starting from a given obj function
840 fFitType = objFunc->Type();
841 fExtObjFunction = nullptr;
842 fObjFunction = std::move(objFunc);
843 if (!DoInitMinimizer()) return false;
844 return DoMinimization(chi2func);
845}
846template<class ObjFunc_t>
847bool Fitter::DoWeightMinimization(std::unique_ptr<ObjFunc_t> objFunc, const ROOT::Math::IMultiGenFunction * chi2func) {
848 // perform the minimization initializing the minimizer starting from a given obj function
849 // and apply afterwards the correction for weights. This applyies only for logL fitting
850 this->fFitType = objFunc->Type();
851 fExtObjFunction = nullptr;
852 fObjFunction = std::move(objFunc);
853 if (!DoInitMinimizer()) return false;
854 if (!DoMinimization(chi2func)) return false;
855 objFunc->UseSumOfWeightSquare();
856 return ApplyWeightCorrection(*objFunc);
857}
858
859
861 // update the fit configuration after a fit using the obtained result
862 if (fResult->IsEmpty() || !fResult->IsValid() ) return;
863 for (unsigned int i = 0; i < fConfig.NPar(); ++i) {
865 par.SetValue( fResult->Value(i) );
866 if (fResult->Error(i) > 0) par.SetStepSize( fResult->Error(i) );
867 }
868}
869
871 // retrieve ncalls from the fit method functions
872 // this function is called when minimizer does not provide a way of returning the nnumber of function calls
873 int ncalls = 0;
874 if (!fUseGradient) {
875 const ROOT::Math::FitMethodFunction * fcn = dynamic_cast<const ROOT::Math::FitMethodFunction *>(fObjFunction.get());
876 if (fcn) ncalls = fcn->NCalls();
877 }
878 else {
880 if (fcn) ncalls = fcn->NCalls();
881 }
882 return ncalls;
883}
884
885
886bool Fitter::ApplyWeightCorrection(const ROOT::Math::IMultiGenFunction & loglw2, bool minimizeW2L) {
887 // apply correction for weight square
888 // Compute Hessian of the loglikelihood function using the sum of the weight squared
889 // This method assumes:
890 // - a fit has been done before and a covariance matrix exists
891 // - the objective function is a likelihood function and Likelihood::UseSumOfWeightSquare()
892 // has been called before
893
894 if (fMinimizer.get() == 0) {
895 MATH_ERROR_MSG("Fitter::ApplyWeightCorrection","Must perform first a fit before applying the correction");
896 return false;
897 }
898
899 unsigned int n = loglw2.NDim();
900 // correct errors for weight squared
901 std::vector<double> cov(n*n);
902 bool ret = fMinimizer->GetCovMatrix(&cov[0] );
903 if (!ret) {
904 MATH_ERROR_MSG("Fitter::ApplyWeightCorrection","Previous fit has no valid Covariance matrix");
905 return false;
906 }
907 // need to use new obj function computed with weight-square
908 std::shared_ptr<ROOT::Math::IMultiGenFunction> objFunc(loglw2.Clone());
909 fObjFunction.swap( objFunc );
910
911 // need to re-initialize the minimizer for the changes applied in the
912 // objective functions
913 if (!DoInitMinimizer()) return false;
914
915 //std::cout << "Running Hesse ..." << std::endl;
916
917 // run eventually before a minimization
918 // ignore its error
919 if (minimizeW2L) fMinimizer->Minimize();
920 // run Hesse on the log-likelihood build using sum of weight squared
921 ret = fMinimizer->Hesse();
922 if (!ret) {
923 MATH_ERROR_MSG("Fitter::ApplyWeightCorrection","Error running Hesse on weight2 likelihood - cannot compute errors");
924 return false;
925 }
926
927 if (fMinimizer->CovMatrixStatus() != 3) {
928 MATH_WARN_MSG("Fitter::ApplyWeightCorrection","Covariance matrix for weighted likelihood is not accurate, the errors may be not reliable");
929 if (fMinimizer->CovMatrixStatus() == 2)
930 MATH_WARN_MSG("Fitter::ApplyWeightCorrection","Covariance matrix for weighted likelihood was forced to be defined positive");
931 if (fMinimizer->CovMatrixStatus() <= 0)
932 // probably should have failed before
933 MATH_ERROR_MSG("Fitter::ApplyWeightCorrection","Covariance matrix for weighted likelihood is not valid !");
934 }
935
936 // get Hessian matrix from weight-square likelihood
937 std::vector<double> hes(n*n);
938 ret = fMinimizer->GetHessianMatrix(&hes[0] );
939 if (!ret) {
940 MATH_ERROR_MSG("Fitter::ApplyWeightCorrection","Error retrieving Hesse on weight2 likelihood - cannot compute errors");
941 return false;
942 }
943
944
945 // perform product of matrix cov * hes * cov
946 // since we do not want to add matrix dependence do product by hand
947 // first do hes * cov
948 std::vector<double> tmp(n*n);
949 for (unsigned int i = 0; i < n; ++i) {
950 for (unsigned int j = 0; j < n; ++j) {
951 for (unsigned int k = 0; k < n; ++k)
952 tmp[i*n+j] += hes[i*n + k] * cov[k*n + j];
953 }
954 }
955 // do multiplication now cov * tmp save result
956 std::vector<double> newCov(n*n);
957 for (unsigned int i = 0; i < n; ++i) {
958 for (unsigned int j = 0; j < n; ++j) {
959 for (unsigned int k = 0; k < n; ++k)
960 newCov[i*n+j] += cov[i*n + k] * tmp[k*n + j];
961 }
962 }
963 // update fit result with new corrected covariance matrix
964 unsigned int k = 0;
965 for (unsigned int i = 0; i < n; ++i) {
966 fResult->fErrors[i] = std::sqrt( newCov[i*(n+1)] );
967 for (unsigned int j = 0; j <= i; ++j)
968 fResult->fCovMatrix[k++] = newCov[i *n + j];
969 }
970
971 // restore previous used objective function
972 fObjFunction.swap( objFunc );
973
974 return true;
975}
976
977
978
980 // return a pointer to the binned data used in the fit
981 // works only for chi2 or binned likelihood fits
982 // thus when the objective function stored is a Chi2Func or a PoissonLikelihood
983 // This also set the model function correctly if it has not been set
984
987
990
991 //MATH_INFO_MSG("Fitter::ExamineFCN","Objective function is not of a known type - FitData and ModelFunction objects are not available");
992 return;
993}
994
995 } // end namespace Fit
996
997} // 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
void SetParamsSettings(unsigned int npar, const double *params, const double *vstep=0)
set the parameter settings from number of parameters and a vector of values and optionally step value...
Definition: FitConfig.cxx:135
bool UpdateAfterFit() const
Update configuration after a fit using the FitResult.
Definition: FitConfig.h:213
void SetMinosErrors(bool on=true)
set Minos erros computation to be performed after fitting
Definition: FitConfig.h:231
void SetMinimizer(const char *type, const char *algo=0)
set minimizer type
Definition: FitConfig.h:181
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
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:551
bool FitFCN()
Perform a fit with the previously set FCN function.
Definition: Fitter.cxx:295
void DoUpdateFitConfig()
Definition: Fitter.cxx:860
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:521
bool DoMinimization(std::unique_ptr< ObjFunc_t > f, const ROOT::Math::IMultiGenFunction *chifunc=nullptr)
do minimization
Definition: Fitter.cxx:838
int fDataSize
size of data sets (need for Fumili or LM fitters)
Definition: Fitter.h:535
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:545
bool DoWeightMinimization(std::unique_ptr< ObjFunc_t > f, const ROOT::Math::IMultiGenFunction *chifunc=nullptr)
Definition: Fitter.cxx:847
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:533
std::shared_ptr< ROOT::Fit::FitData > fData
! pointer to the fit data (binned or unbinned data)
Definition: Fitter.h:547
bool fUseGradient
flag to indicate if using gradient or not
Definition: Fitter.h:527
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:529
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:632
std::shared_ptr< ROOT::Math::IMultiGenFunction > fObjFunction
! pointer to used objective function
Definition: Fitter.h:549
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:886
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:979
const FitConfig & Config() const
access to the fit configuration (const method)
Definition: Fitter.h:410
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:539
std::shared_ptr< ROOT::Fit::FitResult > fResult
! pointer to the object containing the result of the fit
Definition: Fitter.h:543
bool GetDataFromFCN()
internal functions to get data set and model function from FCN useful for fits done with customized F...
Definition: Fitter.h:559
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:537
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:541
bool DoLinearFit()
linear least square fit
Definition: Fitter.cxx:539
bool DoInitMinimizer()
Definition: Fitter.cxx:718
int GetNCallsFromFCN()
Definition: Fitter.cxx:870
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