Logo ROOT  
Reference Guide
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
20#include "TGraph.h" // needed for scan
21#include "TError.h"
22
23#include "TMatrixDSym.h" // needed for inverting the matrix
24
25#include "ThreadLocalStorage.h"
26
27#include <iostream>
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 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 prameter, 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
828 fMinosRun = true;
829
830 double errParab = 0;
831 double gcor = 0;
832 // what returns if parameter fixed or constant or at limit ?
833 fMinuit->mnerrs(i,errUp,errLow, errParab, gcor);
834
835 // do not flag errors case of PROBLEMS or UNCHANGED (
836 return isValid;
837
838}
839
841 // reset TMinuit
842
843 fMinuit->mncler();
844
845 //reset the internal Minuit random generator to its initial state
846 double val = 3;
847 int inseed = 12345;
848 fMinuit->mnrn15(val,inseed);
849
850 fUsed = false;
851 fgUsed = false;
852
853}
854
856 // check if a parameter is defined and in case it was fixed released
857 // TMinuit is not able to release free parameters by redefining them
858 // so we need to force the release
859 if (fMinuit == 0) return;
860 if (fMinuit->GetNumFixedPars() == 0) return;
861 // check if parameter has already been defined
862 if (int(ivar) >= fMinuit->GetNumPars() ) return;
863
864 // check if parameter is fixed
865 for (int i = 0; i < fMinuit->fNpfix; ++i) {
866 if (fMinuit->fIpfix[i] == ivar+1 ) {
867 // parameter is fixed
868 fMinuit->Release(ivar);
869 return;
870 }
871 }
872
873}
874
875
877 // print-out results using classic Minuit format (mnprin)
878 if (fMinuit == 0) return;
879
880 // print minimizer result
881 if (PrintLevel() > 2)
883 else
885}
886
888 // suppress Minuit2 warnings
889 double arglist = 0;
890 int ierr = 0;
891 if (nowarn)
892 fMinuit->mnexcm("SET NOW",&arglist,0,ierr);
893 else
894 fMinuit->mnexcm("SET WAR",&arglist,0,ierr);
895}
896
897
898bool TMinuitMinimizer::Contour(unsigned int ipar, unsigned int jpar, unsigned int &npoints, double * x, double * y) {
899 // contour plot for parameter i and j
900 // need a valid FunctionMinimum otherwise exits
901 if (fMinuit == 0) {
902 Error("TMinuitMinimizer::Contour"," invalid TMinuit instance");
903 return false;
904 }
905
906 // set error and print level
907 double arglist[1];
908 int ierr = 0;
909 arglist[0] = ErrorDef();
910 fMinuit->mnexcm("SET Err",arglist,1,ierr);
911
912 arglist[0] = PrintLevel()-1;
913 fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
914
915 // suppress warning in case Printlevel() == 0
916 if (PrintLevel() == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
917
918 // set precision if needed
919 if (Precision() > 0) {
920 arglist[0] = Precision();
921 fMinuit->mnexcm("SET EPS",arglist,1,ierr);
922 }
923
924
925 if (npoints < 4) {
926 Error("TMinuitMinimizer::Contour","Cannot make contour with so few points");
927 return false;
928 }
929 int npfound = 0;
930 // parameter numbers in mncont start from zero
931 fMinuit->mncont( ipar,jpar,npoints, x, y,npfound);
932 if (npfound<4) {
933 // mncont did go wrong
934 Error("TMinuitMinimizer::Contour","Cannot find more than 4 points");
935 return false;
936 }
937 if (npfound!=(int)npoints) {
938 // mncont did go wrong
939 Warning("TMinuitMinimizer::Contour","Returning only %d points ",npfound);
940 npoints = npfound;
941 }
942 return true;
943
944}
945
946bool TMinuitMinimizer::Scan(unsigned int ipar, unsigned int & nstep, double * x, double * y, double xmin, double xmax) {
947 // scan a parameter (variable) around the minimum value
948 // the parameters must have been set before
949 // if xmin=0 && xmax == 0 by default scan around 2 sigma of the error
950 // if the errors are also zero then scan from min and max of parameter range
951 // (if parameters are limited Minuit scan from min and max instead of 2 sigma by default)
952 // (force in that case to use errors)
953
954 // scan is not implemented for TMinuit, the way to return the array is only via the graph
955 if (fMinuit == 0) {
956 Error("TMinuitMinimizer::Scan"," invalid TMinuit instance");
957 return false;
958 }
959
960 // case of default xmin and xmax
961 if (xmin >= xmax && (int) ipar < fMinuit->GetNumPars() ) {
962 double val = 0; double err = 0;
964 double xlow = 0; double xup = 0 ;
965 int iuint = 0;
966 // in mnpout index starts from ze
967 fMinuit->mnpout( ipar, name, val, err, xlow, xup, iuint);
968 // redefine 2 sigma for all parameters by default (TMinuit does 1 sigma and range if limited)
969 if (iuint > 0 && err > 0) {
970 xmin = val - 2.*err;
971 xmax = val + 2 * err;
972 }
973 }
974
975 double arglist[4];
976 int ierr = 0;
977
978 arglist[0] = PrintLevel()-1;
979 fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
980 // suppress warning in case Printlevel() == 0
981 if (PrintLevel() == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
982
983 // set precision if needed
984 if (Precision() > 0) {
985 arglist[0] = Precision();
986 fMinuit->mnexcm("SET EPS",arglist,1,ierr);
987 }
988
989 if (nstep == 0) return false;
990 arglist[0] = ipar+1; // TMinuit starts from 1
991 arglist[1] = nstep; //+2; // TMinuit deletes two points
992 int nargs = 2;
993 if (xmax > xmin ) {
994 arglist[2] = xmin;
995 arglist[3] = xmax;
996 nargs = 4;
997 }
998 fMinuit->mnexcm("SCAN",arglist,nargs,ierr);
999 if (ierr) {
1000 Error("TMinuitMinimizer::Scan"," Error executing command SCAN");
1001 return false;
1002 }
1003 // get TGraph object
1004 TGraph * gr = dynamic_cast<TGraph *>(fMinuit->GetPlot() );
1005 if (!gr) {
1006 Error("TMinuitMinimizer::Scan"," Error in returned graph object");
1007 return false;
1008 }
1009 nstep = std::min(gr->GetN(), (int) nstep);
1010
1011
1012 std::copy(gr->GetX(), gr->GetX()+nstep, x);
1013 std::copy(gr->GetY(), gr->GetY()+nstep, y);
1014 nstep = gr->GetN();
1015 return true;
1016}
1017
1019 // perform calculation of Hessian
1020
1021 if (fMinuit == 0) {
1022 Error("TMinuitMinimizer::Hesse","invalid TMinuit pointer. Need to call first SetFunction and SetVariable");
1023 return false;
1024 }
1025
1026
1027 double arglist[10];
1028 int ierr = 0;
1029
1030 // set error and print level
1031 arglist[0] = ErrorDef();
1032 fMinuit->mnexcm("SET ERR",arglist,1,ierr);
1033
1034 int printlevel = PrintLevel();
1035 arglist[0] = printlevel - 1;
1036 fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
1037
1038 // suppress warning in case Printlevel() == 0
1039 if (printlevel == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
1040
1041 // set precision if needed
1042 if (Precision() > 0) {
1043 arglist[0] = Precision();
1044 fMinuit->mnexcm("SET EPS",arglist,1,ierr);
1045 }
1046
1047 arglist[0] = MaxFunctionCalls();
1048
1049 fMinuit->mnexcm("HESSE",arglist,1,ierr);
1050 fStatus += 100*ierr;
1051
1052 if (ierr != 0) return false;
1053
1054 // retrieve results (parameter and error matrix)
1055 // only if result is OK
1056
1059
1060 return true;
1061}
1062
1064 // set debug mode
1065
1066 if (fMinuit == 0) {
1067 Error("TMinuitMinimizer::SetDebug","invalid TMinuit pointer. Need to call first SetFunction and SetVariable");
1068 return false;
1069 }
1070 int ierr = 0;
1071 double arglist[1];
1072 arglist[0] = 1;
1073 if (on)
1074 fMinuit->mnexcm("SET DEBUG",arglist,1,ierr);
1075 else
1076 fMinuit->mnexcm("SET NODEBUG",arglist,1,ierr);
1077
1078 return (ierr == 0);
1079}
1080// } // end namespace Fit
1081
1082// } // end namespace ROOT
PyObject * fType
#define f(i)
Definition: RSha256.hxx:104
#define g(i)
Definition: RSha256.hxx:105
#define ClassImp(name)
Definition: Rtypes.h:365
#define R__ASSERT(e)
Definition: TError.h:96
void Info(const char *location, const char *msgfmt,...)
void Error(const char *location, const char *msgfmt,...)
void Warning(const char *location, const char *msgfmt,...)
char name[80]
Definition: TGX11.cxx:109
int type
Definition: TGX11.cxx:120
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
static ROOT::Math::IMultiGenFunction *& GetGlobalFuncPtr()
R__EXTERN TMinuit * gMinuit
Definition: TMinuit.h:271
#define gROOT
Definition: TROOT.h:415
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:327
virtual void Gradient(const T *x, T *grad) const
Evaluate all the vector of function derivatives (gradient) at a point x.
Definition: IFunction.h:342
virtual unsigned int NDim() const=0
Retrieve the dimension of the function.
double Tolerance() const
absolute tolerance
Definition: Minimizer.h:420
unsigned int MaxFunctionCalls() const
max number of function calls
Definition: Minimizer.h:414
double Precision() const
precision of minimizer in the evaluation of the objective function ( a value <=0 corresponds to the l...
Definition: Minimizer.h:424
int Strategy() const
strategy
Definition: Minimizer.h:427
double ErrorDef() const
return the statistical scale used for calculate the error is typically 1 for Chi2 and 0....
Definition: Minimizer.h:434
bool IsValidError() const
return true if Minimizer has performed a detailed error validation (e.g. run Hesse for Minuit)
Definition: Minimizer.h:437
int PrintLevel() const
minimizer configuration parameters
Definition: Minimizer.h:411
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
Double_t * GetY() const
Definition: TGraph.h:131
Int_t GetN() const
Definition: TGraph.h:123
Double_t * GetX() const
Definition: TGraph.h:130
virtual const Element * GetMatrixArray() const
Definition: TMatrixTSym.h:189
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:850
virtual Int_t FixParameter(Int_t parNo)
fix a parameter
Definition: TMinuit.cxx:836
virtual Int_t GetNumPars() const
returns the total number of parameters that have been defined as fixed or free.
Definition: TMinuit.cxx:881
virtual Int_t GetNumFixedPars() const
returns the number of currently fixed parameters
Definition: TMinuit.cxx:864
virtual Int_t Release(Int_t parNo)
release a parameter
Definition: TMinuit.cxx:903
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:1112
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:929
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:1404
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:6256
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:2510
virtual void mnrn15(Double_t &val, Int_t &inseed)
This is a super-portable random number generator.
Definition: TMinuit.cxx:6626
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:2587
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:2673
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:704
virtual void mnprin(Int_t inkode, Double_t fval)
Prints the values of the parameters at the time of the call.
Definition: TMinuit.cxx:6313
Double_t fAmin
Definition: TMinuit.h:49
Double_t fEDM
Definition: TMinuit.h:51
Basic string class.
Definition: TString.h:131
const char * Data() const
Definition: TString.h:364
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
TGraphErrors * gr
Definition: legend1.C:25
RooCmdArg Minimizer(const char *type, const char *alg=0)
auto * m
Definition: textangle.C:8
auto * l
Definition: textangle.C:4