Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TMinuitMinimizer.cxx
Go to the documentation of this file.
1// @(#)root/minuit:$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 TMinuitMinimizer
12
13#include "TMinuitMinimizer.h"
14#include "Math/IFunction.h"
16
17#include "TMinuit.h"
18#include "TROOT.h"
19#include "TList.h"
20
21#include "TGraph.h" // needed for scan
22#include "TError.h"
23
24#include "TMatrixDSym.h" // needed for inverting the matrix
25
26#include "ThreadLocalStorage.h"
27
28#include <cassert>
29#include <algorithm>
30#include <functional>
31#include <cmath>
32
33////////////////////////////////////////////////////////////////////////////////
34/// \class TMinuitMinimizer
35/// \note See ROOT::Minuit2 for a newer version of this class
36/// TMinuitMinimizer class implementing the ROOT::Math::Minimizer interface
37/// using TMinuit. This class is normally instantiated using the plug-in manager
38/// (plug-in with name Minuit or TMinuit).
39/// In addition the user can choose the minimizer algorithm:
40/// Migrad (the default one), Simplex, or Minimize (combined Migrad + Simplex)
41////////////////////////////////////////////////////////////////////////////////
42
43// initialize the static instances
44
45// Implement a thread local static member
48 return fgFunc;
49}
51bool TMinuitMinimizer::fgUsed = false;
52bool TMinuitMinimizer::fgUseStaticMinuit = true; // default case use static Minuit instance
53
54
55
57 fUsed(false),
58 fMinosRun(false),
59 fDim(ndim),
60 fType(type),
61 fMinuit(nullptr)
62{
63 // Constructor for TMinuitMinimier class via an enumeration specifying the minimization
64 // algorithm type. Supported types are : kMigrad, kSimplex, kCombined (a combined
65 // Migrad + Simplex minimization) and kMigradImproved (a Migrad mininimization folloed by an
66 // improved search for global minima). The default type is Migrad (kMigrad).
67
68 // initialize if npar is given
69 if (fDim > 0) InitTMinuit(fDim);
70}
71
72TMinuitMinimizer::TMinuitMinimizer(const char * type, unsigned int ndim ) :
73 fUsed(false),
74 fMinosRun(false),
75 fDim(ndim),
76 fMinuit(nullptr)
77{
78 // constructor from a char * for the algorithm type, used by the plug-in manager
79 // The names supported (case unsensitive) are:
80 // Migrad (default), Simplex, Minimize (for the combined Migrad+ Simplex) and Migrad_imp
81
82 // select type from the string
83 std::string algoname(type);
84 std::transform(algoname.begin(), algoname.end(), algoname.begin(), (int(*)(int)) tolower );
85
87 if (algoname == "simplex") algoType = ROOT::Minuit::kSimplex;
88 if (algoname == "minimize" ) algoType = ROOT::Minuit::kCombined;
89 if (algoname == "migradimproved" ) algoType = ROOT::Minuit::kMigradImproved;
90 if (algoname == "scan" ) algoType = ROOT::Minuit::kScan;
91 if (algoname == "seek" ) algoType = ROOT::Minuit::kSeek;
92
94
95 // initialize if npar is given
96 if (fDim > 0) InitTMinuit(fDim);
97
98}
99
101{
102 // Destructor implementation.
103 if (fMinuit && !fgUseStaticMinuit) {
104 delete fMinuit;
105 fgMinuit = nullptr;
106 }
107}
108
110 // static method to control usage of global TMinuit instance
111 bool prev = fgUseStaticMinuit;
113 return prev;
114}
115
117
118 // when called a second time check dimension - create only if needed
119 // initialize the minuit instance - recreating a new one if needed
120 if (fMinuit ==nullptr || dim > fMinuit->fMaxpar) {
121
122 // case not using the global instance - recreate it all the time
123 if (fgUseStaticMinuit) {
124
125 // re-use gMinuit as static instance of TMinuit
126 // which can be accessed by the user after minimization
127 // check if fgMinuit is different than gMinuit
128 // case 1: fgMinuit not zero but fgMinuit has been deleted (not in gROOT): set to zero
129 // case 2: fgMinuit not zero and exists in global list : set fgMinuit to gMinuit
130 // case 3: fgMinuit zero - and gMinuit not zero: create a new instance locally to avoid conflict
131 if (fgMinuit != gMinuit) {
132 // if object exists in gROOT remove it to avoid a memory leak
133 if (fgMinuit ) {
134 if (gROOT->GetListOfSpecials()->FindObject(fgMinuit) == nullptr) {
135 // case 1: object does not exists in gROOT - means it has been deleted
136 fgMinuit = nullptr;
137 }
138 else {
139 // case 2: object exists - but gMinuit points to something else
140 // restore gMinuit to the one used before by TMinuitMinimizer
142 }
143 }
144 else {
145 // case 3: avoid reusing existing one - maintain fgMinuit to zero
146 // otherwise we will get a double delete if user deletes externally gMinuit
147 // in this case we will loose gMinuit instance
148// fgMinuit = gMinuit;
149// fgUsed = true; // need to reset in case other gMinuit instance is later used
150 }
151 }
152
153 // check if need to create a new TMinuit instance
154 if (fgMinuit == nullptr) {
155 fgUsed = false;
156 fgMinuit = new TMinuit(dim);
157 }
158 else if (fgMinuit->GetNumPars() != int(dim) ) {
159 delete fgMinuit;
160 fgUsed = false;
161 fgMinuit = new TMinuit(dim);
162 }
163
165 }
166
167 else {
168 // re- create all the time a new instance of TMinuit (fgUseStaticMinuit is false)
169 if (fMinuit) delete fMinuit;
170 fMinuit = new TMinuit(dim);
172 fgUsed = false;
173 }
174
175 } // endif fMinuit ==0 || dim > fMaxpar
176
177 fDim = dim;
178
180
181 // set print level in TMinuit
182 double arglist[1];
183 int ierr= 0;
184 // TMinuit level is shift by 1 -1 means 0;
185 arglist[0] = PrintLevel() - 1;
186 fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
188}
189
190
192 // Set the objective function to be minimized, by passing a function object implement the
193 // basic multi-dim Function interface. In this case the derivatives will be
194 // calculated by Minuit
195 // Here a TMinuit instance is created since only at this point we know the number of parameters
196
197 const bool hasGrad = func.HasGradient();
198
199 fDim = func.NDim();
200
201 // create TMinuit if needed
203
204 // assign to the static pointer (NO Thread safety here)
205 GetGlobalFuncPtr() = const_cast<ROOT::Math::IMultiGenFunction *>(&func);
207
208 double arglist[1];
209 int ierr = 0;
210
211 if(hasGrad) {
212 // set gradient
213 // by default do not check gradient calculation
214 // it cannot be done here, check can be done only after having defined the parameters
215 arglist[0] = 1;
216 fMinuit->mnexcm("SET GRAD",arglist,1,ierr);
217 } else {
218 // switch off gradient calculations
219 fMinuit->mnexcm("SET NOGrad",arglist,0,ierr);
220 }
221}
222
223void TMinuitMinimizer::Fcn( int &, double * , double & f, double * x , int /* iflag */) {
224 // implementation of FCN static function used internally by TMinuit.
225 // Adapt IMultiGenFunction interface to TMinuit FCN static function
226 f = GetGlobalFuncPtr()->operator()(x);
227}
228
229void TMinuitMinimizer::FcnGrad( int &, double * g, double & f, double * x , int iflag ) {
230 // implementation of FCN static function used internally by TMinuit.
231 // Adapt IMultiGradFunction interface to TMinuit FCN static function in the case of user
232 // provided gradient.
234
235 assert(gFunc != nullptr);
236 f = gFunc->operator()(x);
237
238 // calculates also derivatives
239 if (iflag == 2) gFunc->Gradient(x,g);
240}
241
242bool TMinuitMinimizer::SetVariable(unsigned int ivar, const std::string & name, double val, double step) {
243 // set a free variable.
244 if (!CheckMinuitInstance()) return false;
245
247
248 // clear after minimization when setting params
249 if (fUsed) DoClear();
250
251 // check if parameter was defined and in case it was fixed, release it
253
254 int iret = fMinuit->DefineParameter(ivar , name.c_str(), val, step, 0., 0. );
255 return (iret == 0);
256}
257
258bool TMinuitMinimizer::SetLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double lower, double upper) {
259 // set a limited variable.
260 if (!CheckMinuitInstance()) return false;
261
263
264 // clear after minimization when setting params
265 if (fUsed) DoClear();
266
267 // check if parameter was defined and in case it was fixed, release it
269
270 int iret = fMinuit->DefineParameter(ivar, name.c_str(), val, step, lower, upper );
271 return (iret == 0);
272}
273
274bool TMinuitMinimizer::SetLowerLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double lower ) {
275 // set a lower limited variable
276 // since is not supported in TMinuit , just use a artificial large value
277 Warning("TMinuitMinimizer::SetLowerLimitedVariable","not implemented - use as upper limit 1.E7 instead of +inf");
278 return SetLimitedVariable(ivar, name, val , step, lower, lower+ 1.E7); // use 1.E7 which will make TMinuit happy
279}
280
281bool TMinuitMinimizer::SetUpperLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double upper ) {
282 // set a upper limited variable
283 // since is not supported in TMinuit , just use a artificial large negative value
284 Warning("TMinuitMinimizer::SetUpperLimitedVariable","not implemented - use as lower limit -1.E7 instead of -inf");
285 return SetLimitedVariable(ivar, name, val , step, upper -1.E7, upper);
286}
287
288
290 // check instance of fMinuit
291 if (fMinuit == nullptr) {
292 Error("TMinuitMinimizer::CheckMinuitInstance","Invalid TMinuit pointer. Need to call first SetFunction");
293 return false;
294 }
295 return true;
296}
297
298bool TMinuitMinimizer::CheckVarIndex(unsigned int ivar) const {
299 // check index of Variable (assume fMinuit exists)
300 if ((int) ivar >= fMinuit->fNu ) {
301 Error("TMinuitMinimizer::CheckVarIndex","Invalid parameter index");
302 return false;
303 }
304 return true;
305}
306
307
308bool TMinuitMinimizer::SetFixedVariable(unsigned int ivar, const std::string & name, double val) {
309 // set a fixed variable.
310 if (!CheckMinuitInstance()) return false;
311
312 // clear after minimization when setting params
314
315 // clear after minimization when setting params
316 if (fUsed) DoClear();
317
318 // put an arbitrary step (0.1*abs(value) otherwise TMinuit consider the parameter as constant
319 // constant parameters are treated differently (they are ignored inside TMinuit and not considered in the
320 // total list of parameters)
321 double step = ( val != 0) ? 0.1 * std::abs(val) : 0.1;
322 int iret = fMinuit->DefineParameter(ivar, name.c_str(), val, step, 0., 0. );
323 if (iret == 0) iret = fMinuit->FixParameter(ivar);
324 return (iret == 0);
325}
326
327bool TMinuitMinimizer::SetVariableValue(unsigned int ivar, double val) {
328 // set the value of an existing variable
329 // parameter must exist or return false
330 if (!CheckMinuitInstance()) return false;
331
332 double arglist[2];
333 int ierr = 0;
334
335 arglist[0] = ivar+1; // TMinuit starts from 1
336 arglist[1] = val;
337 fMinuit->mnexcm("SET PAR",arglist,2,ierr);
338 return (ierr==0);
339}
340
341bool TMinuitMinimizer::SetVariableStepSize(unsigned int ivar, double step) {
342 // set the step-size of an existing variable
343 // parameter must exist or return false
344 if (!CheckMinuitInstance()) return false;
345 // need to re-implement re-calling mnparm
346 // get first current parameter values and limits
347 double curval,err, lowlim, uplim;
348 int iuint; // internal index
351 if (iuint == -1) return false;
353 return (iret == 0);
354
355}
356
358 // set the limits of an existing variable
359 // parameter must exist or return false
360 Warning("TMinuitMinimizer::SetVariableLowerLimit","not implemented - use as upper limit 1.E30 instead of +inf");
361 return SetVariableLimits(ivar, lower, 1.E30);
362}
364 // set the limits of an existing variable
365 // parameter must exist or return false
366 Warning("TMinuitMinimizer::SetVariableUpperLimit","not implemented - - use as lower limit -1.E30 instead of +inf");
367 return SetVariableLimits(ivar, -1.E30, upper);
368}
369
370bool TMinuitMinimizer::SetVariableLimits(unsigned int ivar, double lower, double upper) {
371 // set the limits of an existing variable
372 // parameter must exist or return false
373
374 if (!CheckMinuitInstance()) return false;
375 // need to re-implement re-calling mnparm
376 // get first current parameter values and limits
377 double curval,err, lowlim, uplim;
378 int iuint; // internal index
381 if (iuint == -1) return false;
383 return (iret == 0);
384
385}
386
388 // Fix an existing variable
389 if (!CheckMinuitInstance()) return false;
390 if (!CheckVarIndex(ivar)) return false;
392 return (iret == 0);
393}
394
396 // Fix an existing variable
397 if (!CheckMinuitInstance()) return false;
398 if (!CheckVarIndex(ivar)) return false;
399 int iret = fMinuit->Release(ivar);
400 return (iret == 0);
401}
402
403bool TMinuitMinimizer::IsFixedVariable(unsigned int ivar) const {
404 // query if variable is fixed
405 if (!CheckMinuitInstance()) return false;
406 if (!CheckVarIndex(ivar)) return false;
407 return (fMinuit->fNiofex[ivar] == 0 );
408}
409
411 // retrieve variable settings (all set info on the variable)
412 if (!CheckMinuitInstance()) return false;
413 if (!CheckVarIndex(ivar)) return false;
414 double curval,err, lowlim, uplim;
415 int iuint; // internal index
418 if (iuint == -1) return false;
419 var.Set(name.Data(), curval, err, lowlim, uplim);
420 if (IsFixedVariable(ivar)) var.Fix();
421 return true;
422}
423
424
425
426std::string TMinuitMinimizer::VariableName(unsigned int ivar) const {
427 // return the variable name
428 if (!CheckMinuitInstance()) return std::string();
429 if (!CheckVarIndex(ivar)) return std::string();
430 return fMinuit->fCpnam[ivar].Data();
431}
432
433int TMinuitMinimizer::VariableIndex(const std::string & ) const {
434 // return variable index
435 Error("TMinuitMinimizer::VariableIndex"," find index of a variable from its name is not implemented in TMinuit");
436 return -1;
437}
438
440 // perform the minimization using the algorithm chosen previously by the user
441 // By default Migrad is used.
442 // Return true if the found minimum is valid and update internal cached values of
443 // minimum values, errors and covariance matrix.
444 // Status of minimizer is set to:
445 // migradResult + 10*minosResult + 100*hesseResult + 1000*improveResult
446
447
448 if (fMinuit == nullptr) {
449 Error("TMinuitMinimizer::Minimize","invalid TMinuit pointer. Need to call first SetFunction and SetVariable");
450 return false;
451 }
452
453
454 // total number of parameter defined in Minuit is fNu
455 if (fMinuit->fNu < static_cast<int>(fDim) ) {
456 Error("TMinuitMinimizer::Minimize","The total number of defined parameters is different than the function dimension, npar = %d, dim = %d",fMinuit->fNu, fDim);
457 return false;
458 }
459
460 int printlevel = PrintLevel();
461
462 // total number of free parameter is 0
463 if (fMinuit->fNpar <= 0) {
464 // retrieve parameters values from TMinuit
466 fMinuit->fAmin = (*GetGlobalFuncPtr())(&fParams.front());
467 if (printlevel > 0) Info("TMinuitMinimizer::Minimize","There are no free parameter - just compute the function value");
468 return true;
469 }
470
471
472 double arglist[10];
473 int ierr = 0;
474
475
476 // set error and print level
477 arglist[0] = ErrorDef();
478 fMinuit->mnexcm("SET Err",arglist,1,ierr);
479
480 arglist[0] = printlevel - 1;
481 fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
482
483 // suppress warning in case Printlevel() == 0
484 if (printlevel == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
485
486 // set precision if needed
487 if (Precision() > 0) {
488 arglist[0] = Precision();
489 fMinuit->mnexcm("SET EPS",arglist,1,ierr);
490 }
491
492 // set strategy
493 int strategy = Strategy();
494 if (strategy >=0 && strategy <=2 ) {
495 arglist[0] = strategy;
496 fMinuit->mnexcm("SET STR",arglist,1,ierr);
497 }
498
500 arglist[1] = Tolerance();
501
502 int nargs = 2;
503
504 switch (fType){
506 // case of Migrad
507 fMinuit->mnexcm("MIGRAD",arglist,nargs,ierr);
508 break;
510 // case of combined (Migrad+ simplex)
511 fMinuit->mnexcm("MINIMIZE",arglist,nargs,ierr);
512 break;
514 // case of Simlex
515 fMinuit->mnexcm("SIMPLEX",arglist,nargs,ierr);
516 break;
518 // case of Scan (scan all parameters with default values)
519 nargs = 0;
521 break;
523 // case of Seek (random find minimum in a hypercube around current parameter values
524 // use Tolerance as measures for standard deviation (if < 1) used default value in Minuit ( supposed to be 3)
525 nargs = 1;
526 if (arglist[1] >= 1.) nargs = 2;
528 break;
529 default:
530 // default: use Migrad
531 fMinuit->mnexcm("MIGRAD",arglist,nargs,ierr);
532
533 }
534
535 fgUsed = true;
536 fUsed = true;
537
538 fStatus = ierr;
539 int minErrStatus = ierr;
540
541 if (printlevel>2) Info("TMinuitMinimizer::Minimize","Finished to run MIGRAD - status %d",ierr);
542
543 // run improved if needed
545 fMinuit->mnexcm("IMPROVE",arglist,1,ierr);
546 fStatus += 1000*ierr;
547 if (printlevel>2) Info("TMinuitMinimizer::Minimize","Finished to run IMPROVE - status %d",ierr);
548 }
549
550
551 // check if Hesse needs to be run
552 // Migrad runs inside it automatically for strategy >=1. Do also
553 // in case improve or other minimizers are used
554 // run Hesse also in case cov matrix has not been computed or has been made pos-def
555 // or a valid error analysis is requested (when IsValidError() == true)
556 if (minErrStatus == 0 && (IsValidError() || ( strategy >=1 && CovMatrixStatus() < 3) ) ) {
557 fMinuit->mnexcm("HESSE",arglist,1,ierr);
558 fStatus += 100*ierr;
559 // here ierr is zero when Hesse fails CovMatrixStatus = 0 or 1 (2 is when was made posdef)
560 if (ierr == 0 && CovMatrixStatus() < 2){
561 fStatus += 100*(CovMatrixStatus()+1);
562 }
563 // should check here cov matrix status??? (if <3 flag error ?)
564 if (printlevel>2) Info("TMinuitMinimizer::Minimize","Finished to run HESSE - status %d",ierr);
565 }
566
567 // retrieve parameters and errors from TMinuit
569
570 if (minErrStatus == 0) {
571
572 // store global min results (only if minimization is OK)
573 // ignore cases when Hesse or IMprove return error different than zero
575
576 // need to re-run Minos again if requested
577 fMinosRun = false;
578
579 return true;
580
581 }
582 return false;
583
584}
585
587 // retrieve from TMinuit minimum parameter values
588 // and errors
589
590 assert(fMinuit != nullptr);
591
592 // get parameter values
593 if (fParams.size() != fDim) fParams.resize( fDim);
594 if (fErrors.size() != fDim) fErrors.resize( fDim);
595 for (unsigned int i = 0; i < fDim; ++i) {
596 fMinuit->GetParameter( i, fParams[i], fErrors[i]);
597 }
598}
599
601 // get covariance error matrix from TMinuit
602 // when some parameters are fixed filled the corresponding rows and column with zero's
603
604 assert(fMinuit != nullptr);
605
606 unsigned int nfree = NFree();
607
608 unsigned int ndim2 = fDim*fDim;
609 if (fCovar.size() != ndim2 ) fCovar.resize(fDim*fDim);
610 if (nfree >= fDim) { // no fixed parameters
611 fMinuit->mnemat(&fCovar.front(), fDim);
612 }
613 else {
614 // case of fixed params need to take care
615 std::vector<double> tmpMat(nfree*nfree);
616 fMinuit->mnemat(&tmpMat.front(), nfree);
617
618 unsigned int l = 0;
619 for (unsigned int i = 0; i < fDim; ++i) {
620
621 if ( fMinuit->fNiofex[i] > 0 ) { // not fixed ?
622 unsigned int m = 0;
623 for (unsigned int j = 0; j <= i; ++j) {
624 if ( fMinuit->fNiofex[j] > 0 ) { //not fixed
625 fCovar[i*fDim + j] = tmpMat[l*nfree + m];
626 fCovar[j*fDim + i] = fCovar[i*fDim + j];
627 m++;
628 }
629 }
630 l++;
631 }
632 }
633
634 }
635}
636
637unsigned int TMinuitMinimizer::NCalls() const {
638 // return total number of function calls
639 if (fMinuit == nullptr) return 0;
640 return fMinuit->fNfcn;
641}
642
644 // return minimum function value
645
646 // use part of code from mnstat
647 if (!fMinuit) return 0;
648 double minval = fMinuit->fAmin;
649 if (minval == fMinuit->fUndefi) return 0;
650 return minval;
651}
652
653double TMinuitMinimizer::Edm() const {
654 // return expected distance from the minimum
655
656 // use part of code from mnstat
657 if (!fMinuit) return -1;
658 if (fMinuit->fAmin == fMinuit->fUndefi || fMinuit->fEDM == fMinuit->fBigedm) return fMinuit->fUp;
659 return fMinuit->fEDM;
660}
661
662unsigned int TMinuitMinimizer::NFree() const {
663 // return number of free parameters
664 if (!fMinuit) return 0;
665 if (fMinuit->fNpar < 0) return 0;
666 return fMinuit->fNpar;
667}
668
670 // get covariance matrix
672 if ( fCovar.size() != fDim*fDim || covStatus < 2) {
673 Error("TMinuitMinimizer::GetHessianMatrix","Hessian matrix has not been computed - status %d",covStatus);
674 return false;
675 }
676 std::copy(fCovar.begin(), fCovar.end(), cov);
678 return true;
679}
680
682 // get Hessian - inverse of covariance matrix
683 // just invert it
684 // but need to get the compact form to avoid the zero for the fixed parameters
686 if ( fCovar.size() != fDim*fDim || covStatus < 2) {
687 Error("TMinuitMinimizer::GetHessianMatrix","Hessian matrix has not been computed - status %d",covStatus);
688 return false;
689 }
690 // case of fixed params need to take care
691 unsigned int nfree = NFree();
693 fMinuit->mnemat(mat.GetMatrixArray(), nfree);
694 // invert the matrix
695 // probably need to check if failed. In that case inverse is equal to original
696 mat.Invert();
697
698 unsigned int l = 0;
699 for (unsigned int i = 0; i < fDim; ++i) {
700 if ( fMinuit->fNiofex[i] > 0 ) { // not fixed ?
701 unsigned int m = 0;
702 for (unsigned int j = 0; j <= i; ++j) {
703 if ( fMinuit->fNiofex[j] > 0 ) { //not fixed
704 hes[i*fDim + j] = mat(l,m);
705 hes[j*fDim + i] = hes[i*fDim + j];
706 m++;
707 }
708 }
709 l++;
710 }
711 }
712 return true;
713}
714// if ( fCovar.size() != fDim*fDim ) return false;
715// TMatrixDSym mat(fDim, &fCovar.front() );
716// std::copy(mat.GetMatrixArray(), mat.GetMatrixArray()+ mat.GetNoElements(), hes);
717// return true;
718// }
719
721 // return status of covariance matrix
722 // status: 0= not calculated at all
723 // 1= approximation only, not accurate
724 // 2= full matrix, but forced positive-definite
725 // 3= full accurate covariance matrix
726
727 // use part of code from mnstat
728 if (!fMinuit) return 0;
729 if (fMinuit->fAmin == fMinuit->fUndefi) return 0;
730 return fMinuit->fISW[1];
731}
732
733double TMinuitMinimizer::GlobalCC(unsigned int i) const {
734 // global correlation coefficient for parameter i
735 if (!fMinuit) return 0;
736 if (!fMinuit->fGlobcc) return 0;
737 if (int(i) >= fMinuit->fNu) return 0;
738 // get internal number in Minuit
739 int iin = fMinuit->fNiofex[i];
740 // index in TMinuit starts from 1
741 if (iin < 1) return 0;
742 return fMinuit->fGlobcc[iin-1];
743}
744
745bool TMinuitMinimizer::GetMinosError(unsigned int i, double & errLow, double & errUp, int ) {
746 // Perform Minos analysis for the given parameter i
747
748 if (fMinuit == nullptr) {
749 Error("TMinuitMinimizer::GetMinosError","invalid TMinuit pointer. Need to call first SetFunction and SetVariable");
750 return false;
751 }
752
753 // check if parameter is fixed
754 if (fMinuit->fNiofex[i] == 0 ) {
755 if (PrintLevel() > 0) Info("TMinuitMinimizer::GetMinosError","Parameter %s is fixed. There are no Minos error to calculate. Ignored.",VariableName(i).c_str());
756 errLow = 0; errUp = 0;
757 return true;
758 }
759
760 double arglist[2];
761 int ierr = 0;
762
763
764 // set error, print level, precision and strategy if they have changed
765 if (fMinuit->fUp != ErrorDef() ) {
766 arglist[0] = ErrorDef();
767 fMinuit->mnexcm("SET Err",arglist,1,ierr);
768 }
769
770 if (fMinuit->fISW[4] != (PrintLevel()-1) ) {
771 arglist[0] = PrintLevel()-1;
772 fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
773 // suppress warning in case Printlevel() == 0
774 if (PrintLevel() == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
775 }
776 if (fMinuit->fIstrat != Strategy() ) {
777 arglist[0] = Strategy();
778 fMinuit->mnexcm("SET STR",arglist,1,ierr);
779 }
780
781 if (Precision() > 0 && fMinuit->fEpsma2 != Precision() ) {
782 arglist[0] = Precision();
783 fMinuit->mnexcm("SET EPS",arglist,1,ierr);
784 }
785
786
787 // syntax of MINOS command is: MINOS [maxcalls] [parno]
788 // if parno = 0 all parameters are done
790 arglist[1] = i+1; // par number starts from 1 in TMInuit
791
792 int nargs = 2;
793 fMinuit->mnexcm("MINOS",arglist,nargs,ierr);
794 bool isValid = (ierr == 0);
795 // check also the status from fCstatu
796 if (isValid && fMinuit->fCstatu != "SUCCESSFUL") {
797 if (fMinuit->fCstatu == "FAILURE" ) {
798 // in this case MINOS failed on all parameters, so it is not valid !
799 ierr = 5;
800 isValid = false;
801 }
802 if (fMinuit->fCstatu == "PROBLEMS") ierr = 6;
803 ierr = 7; // this should be the case UNCHANGED
804 }
805
806 fStatus += 10*ierr;
808
809 fMinosRun = true;
810
811 // retrieve parameters in case a new minimum has been found
812 if (fMinuit->fCstatu == "SUCCESSFUL")
814
815 double errParab = 0;
816 double gcor = 0;
817 // what returns if parameter fixed or constant or at limit ?
819
820 // do not flag errors case of PROBLEMS or UNCHANGED (
821 return isValid;
822
823}
824
826 // reset TMinuit
827
828 fMinuit->mncler();
829
830 //reset the internal Minuit random generator to its initial state
831 double val = 3;
832 int inseed = 12345;
833 fMinuit->mnrn15(val,inseed);
834
835 fUsed = false;
836 fgUsed = false;
837
838}
839
841 // check if a parameter is defined and in case it was fixed released
842 // TMinuit is not able to release free parameters by redefining them
843 // so we need to force the release
844 if (fMinuit == nullptr) return;
845 if (fMinuit->GetNumFixedPars() == 0) return;
846 // check if parameter has already been defined
847 if (int(ivar) >= fMinuit->GetNumPars() ) return;
848
849 // check if parameter is fixed
850 for (int i = 0; i < fMinuit->fNpfix; ++i) {
851 if (fMinuit->fIpfix[i] == ivar+1 ) {
852 // parameter is fixed
854 return;
855 }
856 }
857
858}
859
860
862 // print-out results using classic Minuit format (mnprin)
863 if (fMinuit == nullptr) return;
864
865 // print minimizer result
866 if (PrintLevel() > 2)
868 else
870}
871
873 // suppress Minuit2 warnings
874 double arglist = 0;
875 int ierr = 0;
876 if (nowarn)
877 fMinuit->mnexcm("SET NOW",&arglist,0,ierr);
878 else
879 fMinuit->mnexcm("SET WAR",&arglist,0,ierr);
880}
881
882
883bool TMinuitMinimizer::Contour(unsigned int ipar, unsigned int jpar, unsigned int &npoints, double * x, double * y) {
884 // contour plot for parameter i and j
885 // need a valid FunctionMinimum otherwise exits
886 if (fMinuit == nullptr) {
887 Error("TMinuitMinimizer::Contour"," invalid TMinuit instance");
888 return false;
889 }
890
891 // set error and print level
892 double arglist[1];
893 int ierr = 0;
894 arglist[0] = ErrorDef();
895 fMinuit->mnexcm("SET Err",arglist,1,ierr);
896
897 arglist[0] = PrintLevel()-1;
898 fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
899
900 // suppress warning in case Printlevel() == 0
901 if (PrintLevel() == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
902
903 // set precision if needed
904 if (Precision() > 0) {
905 arglist[0] = Precision();
906 fMinuit->mnexcm("SET EPS",arglist,1,ierr);
907 }
908
909
910 if (npoints < 4) {
911 Error("TMinuitMinimizer::Contour","Cannot make contour with so few points");
912 return false;
913 }
914 int npfound = 0;
915 // parameter numbers in mncont start from zero
916 fMinuit->mncont( ipar,jpar,npoints, x, y,npfound);
917 if (npfound<4) {
918 // mncont did go wrong
919 Error("TMinuitMinimizer::Contour","Cannot find more than 4 points");
920 return false;
921 }
922 if (npfound!=(int)npoints) {
923 // mncont did go wrong
924 Warning("TMinuitMinimizer::Contour","Returning only %d points ",npfound);
926 }
927 return true;
928
929}
930
931bool TMinuitMinimizer::Scan(unsigned int ipar, unsigned int & nstep, double * x, double * y, double xmin, double xmax) {
932 // scan a parameter (variable) around the minimum value
933 // the parameters must have been set before
934 // if xmin=0 && xmax == 0 by default scan around 2 sigma of the error
935 // if the errors are also zero then scan from min and max of parameter range
936 // (if parameters are limited Minuit scan from min and max instead of 2 sigma by default)
937 // (force in that case to use errors)
938
939 // scan is not implemented for TMinuit, the way to return the array is only via the graph
940 if (fMinuit == nullptr) {
941 Error("TMinuitMinimizer::Scan"," invalid TMinuit instance");
942 return false;
943 }
944
945 // case of default xmin and xmax
946 if (xmin >= xmax && (int) ipar < fMinuit->GetNumPars() ) {
947 double val = 0; double err = 0;
949 double xlow = 0; double xup = 0 ;
950 int iuint = 0;
951 // in mnpout index starts from ze
952 fMinuit->mnpout( ipar, name, val, err, xlow, xup, iuint);
953 // redefine 2 sigma for all parameters by default (TMinuit does 1 sigma and range if limited)
954 if (iuint > 0 && err > 0) {
955 xmin = val - 2.*err;
956 xmax = val + 2 * err;
957 }
958 }
959
960 double arglist[4];
961 int ierr = 0;
962
963 arglist[0] = PrintLevel()-1;
964 fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
965 // suppress warning in case Printlevel() == 0
966 if (PrintLevel() == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
967
968 // set precision if needed
969 if (Precision() > 0) {
970 arglist[0] = Precision();
971 fMinuit->mnexcm("SET EPS",arglist,1,ierr);
972 }
973
974 if (nstep == 0) return false;
975 arglist[0] = ipar+1; // TMinuit starts from 1
976 arglist[1] = nstep; //+2; // TMinuit deletes two points
977 int nargs = 2;
978 if (xmax > xmin ) {
979 arglist[2] = xmin;
980 arglist[3] = xmax;
981 nargs = 4;
982 }
984 if (ierr) {
985 Error("TMinuitMinimizer::Scan"," Error executing command SCAN");
986 return false;
987 }
988 // get TGraph object
989 TGraph * gr = dynamic_cast<TGraph *>(fMinuit->GetPlot() );
990 if (!gr) {
991 Error("TMinuitMinimizer::Scan"," Error in returned graph object");
992 return false;
993 }
994 nstep = std::min(gr->GetN(), (int) nstep);
995
996
997 std::copy(gr->GetX(), gr->GetX()+nstep, x);
998 std::copy(gr->GetY(), gr->GetY()+nstep, y);
999 nstep = gr->GetN();
1000 return true;
1001}
1002
1004 // perform calculation of Hessian
1005
1006 if (fMinuit == nullptr) {
1007 Error("TMinuitMinimizer::Hesse","invalid TMinuit pointer. Need to call first SetFunction and SetVariable");
1008 return false;
1009 }
1010
1011
1012 double arglist[10];
1013 int ierr = 0;
1014
1015 // set error and print level
1016 arglist[0] = ErrorDef();
1017 fMinuit->mnexcm("SET ERR",arglist,1,ierr);
1018
1019 int printlevel = PrintLevel();
1020 arglist[0] = printlevel - 1;
1021 fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
1022
1023 // suppress warning in case Printlevel() == 0
1024 if (printlevel == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
1025
1026 // set precision if needed
1027 if (Precision() > 0) {
1028 arglist[0] = Precision();
1029 fMinuit->mnexcm("SET EPS",arglist,1,ierr);
1030 }
1031
1033
1034 fMinuit->mnexcm("HESSE",arglist,1,ierr);
1035 fStatus += 100*ierr;
1036
1037 if (ierr != 0) return false;
1038
1039 // retrieve results (parameter and error matrix)
1040 // only if result is OK
1041
1044
1045 return true;
1046}
1047
1049 // set debug mode
1050
1051 if (fMinuit == nullptr) {
1052 Error("TMinuitMinimizer::SetDebug","invalid TMinuit pointer. Need to call first SetFunction and SetVariable");
1053 return false;
1054 }
1055 int ierr = 0;
1056 double arglist[1];
1057 arglist[0] = 1;
1058 if (on)
1059 fMinuit->mnexcm("SET DEBUG",arglist,1,ierr);
1060 else
1061 fMinuit->mnexcm("SET NODEBUG",arglist,1,ierr);
1062
1063 return (ierr == 0);
1064}
1065// } // end namespace Fit
1066
1067// } // end namespace ROOT
#define f(i)
Definition RSha256.hxx:104
#define g(i)
Definition RSha256.hxx:105
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition TError.cxx:241
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:208
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition TError.cxx:252
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void on
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
static ROOT::Math::IMultiGenFunction *& GetGlobalFuncPtr()
R__EXTERN TMinuit * gMinuit
Definition TMinuit.h:271
#define gROOT
Definition TROOT.h:411
Class, describing value, limits and step size of the parameters Provides functionality also to set/re...
void Set(const std::string &name, double value, double step)
set value and name (unlimited parameter)
void Fix()
fix the parameter
Documentation for the abstract class IBaseFunctionMultiDim.
Definition IFunction.h:63
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:252
double Tolerance() const
Absolute tolerance.
Definition Minimizer.h:317
unsigned int MaxFunctionCalls() const
Max number of function calls.
Definition Minimizer.h:311
double Precision() const
Precision of minimizer in the evaluation of the objective function.
Definition Minimizer.h:321
int fStatus
status of minimizer
Definition Minimizer.h:388
int Strategy() const
Strategy.
Definition Minimizer.h:324
double ErrorDef() const
Definition Minimizer.h:334
bool IsValidError() const
Definition Minimizer.h:337
int PrintLevel() const
Set print level.
Definition Minimizer.h:308
const_iterator begin() const
const_iterator end() const
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
Double_t * GetY() const
Definition TGraph.h:139
Int_t GetN() const
Definition TGraph.h:131
Double_t * GetX() const
Definition TGraph.h:138
bool FixVariable(unsigned int) override
fix an existing variable
TMinuitMinimizer(ROOT::Minuit::EMinimizerType type=ROOT::Minuit::kMigrad, unsigned int ndim=0)
Default constructor.
void RetrieveErrorMatrix()
retrieve error matrix from TMinuit
bool SetFixedVariable(unsigned int, const std::string &, double) override
set fixed variable (override if minimizer supports them )
static TMinuit * fgMinuit
bool SetVariableLimits(unsigned int ivar, double lower, double upper) override
set the limits of an existing variable
static void Fcn(int &, double *, double &f, double *, int)
implementation of FCN for Minuit
bool SetVariable(unsigned int ivar, const std::string &name, double val, double step) override
set free variable
bool ReleaseVariable(unsigned int) override
release an existing variable
static void FcnGrad(int &, double *g, double &f, double *, int)
implementation of FCN for Minuit when user provided gradient is used
ROOT::Minuit::EMinimizerType fType
bool SetVariableLowerLimit(unsigned int, double) override
set the lower-limit of an existing variable
int VariableIndex(const std::string &name) const override
get index of variable given a variable given a name return always -1 .
double MinValue() const override
return minimum function value
bool SetVariableStepSize(unsigned int, double) override
set the step size of an existing variable
void RetrieveParams()
retrieve minimum parameters and errors from TMinuit
~TMinuitMinimizer() override
Destructor (no operations)
std::vector< double > fErrors
void SetFunction(const ROOT::Math::IMultiGenFunction &func) override
set the function to minimize
bool CheckMinuitInstance() const
check TMinuit instance
static bool fgUseStaticMinuit
std::vector< double > fCovar
bool Hesse() override
perform a full calculation of the Hessian matrix for error calculation
bool CheckVarIndex(unsigned int ivar) const
check parameter
bool Minimize() override
method to perform the minimization
bool SetVariableValue(unsigned int, double) override
set the value of an existing variable
int CovMatrixStatus() const override
return status of covariance matrix
void SuppressMinuitWarnings(bool nowarn=true)
suppress the minuit warnings (if called with false will enable them) By default they are suppressed o...
bool SetVariableUpperLimit(unsigned int, double) override
set the upper-limit of an existing variable
double Edm() const override
return expected distance reached from the minimum
static bool UseStaticMinuit(bool on=true)
static function to switch on/off usage of static global TMinuit instance (gMinuit) By default it is u...
void PrintResults() override
Print the result according to set level (implemented for TMinuit for maintaining Minuit-style printin...
unsigned int NCalls() const override
number of function calls to reach the minimum
void DoReleaseFixParameter(int ivar)
release a parameter that is fixed when it is redefined
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 )
bool GetMinosError(unsigned int i, double &errLow, double &errUp, int=0) override
minos error for variable i, return false if Minos failed
bool GetVariableSettings(unsigned int, ROOT::Fit::ParameterSettings &) const override
get variable settings in a variable object (like ROOT::Fit::ParamsSettings)
double GlobalCC(unsigned int) const override
global correlation coefficient for variable i
bool IsFixedVariable(unsigned int) const override
query if an existing variable is fixed (i.e.
std::string VariableName(unsigned int ivar) const override
return reference to the objective function virtual const ROOT::Math::IGenFunction & Function() const;
bool SetUpperLimitedVariable(unsigned int ivar, const std::string &name, double val, double step, double upper) override
set upper limit variable (override if minimizer supports them )
bool SetLowerLimitedVariable(unsigned int ivar, const std::string &name, double val, double step, double lower) override
set lower limit variable (override if minimizer supports them )
unsigned int NFree() const override
number of free variables (real dimension of the problem) this is <= Function().NDim() which is the to...
std::vector< double > fParams
void InitTMinuit(int ndim)
initialize the TMinuit instance
bool GetHessianMatrix(double *h) const override
Fill the passed array with the Hessian matrix elements The Hessian matrix is the matrix of the second...
bool Scan(unsigned int i, unsigned int &nstep, double *x, double *y, double xmin=0, double xmax=0) override
scan a parameter i around the minimum.
bool SetDebug(bool on=true)
set debug mode. Return true if setting was successful
bool GetCovMatrix(double *cov) const override
Fill the passed array with the covariance matrix elements if the variable is fixed or const the value...
bool Contour(unsigned int i, unsigned int j, unsigned int &npoints, double *xi, double *xj) override
find the contour points (xi,xj) of the function for parameter i and j around the minimum The contour ...
virtual Int_t GetParameter(Int_t parNo, Double_t &currentValue, Double_t &currentError) const
return parameter value and error
Definition TMinuit.cxx:840
virtual Int_t FixParameter(Int_t parNo)
fix a parameter
Definition TMinuit.cxx:826
virtual Int_t GetNumPars() const
returns the total number of parameters that have been defined as fixed or free.
Definition TMinuit.cxx:871
virtual Int_t GetNumFixedPars() const
returns the number of currently fixed parameters
Definition TMinuit.cxx:854
virtual Int_t Release(Int_t parNo)
release a parameter
Definition TMinuit.cxx:893
Double_t fUndefi
Definition TMinuit.h:60
Int_t fNu
Definition TMinuit.h:130
Double_t * fGlobcc
Definition TMinuit.h:74
virtual void mncler()
Resets the parameter list to UNDEFINED.
Definition TMinuit.cxx:1102
Double_t fUp
Definition TMinuit.h:50
TString * fCpnam
Character to be plotted at the X,Y contour positions.
Definition TMinuit.h:165
Int_t fNfcn
Definition TMinuit.h:145
Double_t fBigedm
Definition TMinuit.h:61
TString fCstatu
Definition TMinuit.h:167
virtual TObject * GetPlot() const
Definition TMinuit.h:200
Int_t fNpfix
Definition TMinuit.h:37
Int_t fISW[7]
Definition TMinuit.h:141
Int_t * fIpfix
Definition TMinuit.h:129
virtual void SetFCN(void(*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
To set the address of the minimization function.
Definition TMinuit.cxx:919
virtual void mncont(Int_t ke1, Int_t ke2, Int_t nptu, Double_t *xptu, Double_t *yptu, Int_t &ierrf)
Find points along a contour where FCN is minimum.
Definition TMinuit.cxx:1394
Double_t fEpsma2
Definition TMinuit.h:57
virtual void mnpout(Int_t iuext, TString &chnam, Double_t &val, Double_t &err, Double_t &xlolim, Double_t &xuplim, Int_t &iuint) const
Provides the user with information concerning the current status.
Definition TMinuit.cxx:6241
Int_t fIstrat
Definition TMinuit.h:150
virtual void mnemat(Double_t *emat, Int_t ndim)
Calculates the external error matrix from the internal matrix.
Definition TMinuit.cxx:2500
virtual void mnrn15(Double_t &val, Int_t &inseed)
This is a super-portable random number generator.
Definition TMinuit.cxx:6613
virtual void mnerrs(Int_t number, Double_t &eplus, Double_t &eminus, Double_t &eparab, Double_t &gcc)
Utility routine to get MINOS errors.
Definition TMinuit.cxx:2577
Int_t fNpar
Definition TMinuit.h:41
virtual void mnexcm(const char *comand, Double_t *plist, Int_t llist, Int_t &ierflg)
Interprets a command and takes appropriate action.
Definition TMinuit.cxx:2663
Int_t fMaxpar
Definition TMinuit.h:39
Int_t * fNiofex
Definition TMinuit.h:127
virtual Int_t DefineParameter(Int_t parNo, const char *name, Double_t initVal, Double_t initErr, Double_t lowerLimit, Double_t upperLimit)
Define a parameter.
Definition TMinuit.cxx:694
virtual void mnprin(Int_t inkode, Double_t fval)
Prints the values of the parameters at the time of the call.
Definition TMinuit.cxx:6298
Double_t fAmin
Definition TMinuit.h:49
Double_t fEDM
Definition TMinuit.h:51
Basic string class.
Definition TString.h:138
const char * Data() const
Definition TString.h:384
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
TGraphErrors * gr
Definition legend1.C:25
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4