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