Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
RooMinimizerFcn.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitCore *
4  * @(#)root/roofitcore:$Id$
5  * Authors: *
6  * AL, Alfio Lazzaro, INFN Milan, alfio.lazzaro@mi.infn.it *
7  * *
8  * *
9  * Redistribution and use in source and binary forms, *
10  * with or without modification, are permitted according to the terms *
11  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
12  *****************************************************************************/
13 
14 #ifndef __ROOFIT_NOROOMINIMIZER
15 
16 //////////////////////////////////////////////////////////////////////////////
17 //
18 // RooMinimizerFcn is am interface class to the ROOT::Math function
19 // for minization.
20 //
21 
22 #include <iostream>
23 
24 #include "RooFit.h"
25 #include "RooMinimizerFcn.h"
26 
27 #include "Riostream.h"
28 
29 #include "TIterator.h"
30 #include "TClass.h"
31 
32 #include "RooAbsArg.h"
33 #include "RooAbsPdf.h"
34 #include "RooArgSet.h"
35 #include "RooRealVar.h"
36 #include "RooAbsRealLValue.h"
37 #include "RooMsgService.h"
38 
39 #include "RooMinimizer.h"
40 
41 using namespace std;
42 
43 RooMinimizerFcn::RooMinimizerFcn(RooAbsReal *funct, RooMinimizer* context,
44  bool verbose) :
45  _funct(funct), _context(context),
46  // Reset the *largest* negative log-likelihood value we have seen so far
47  _maxFCN(-1e30), _numBadNLL(0),
48  _printEvalErrors(10), _doEvalErrorWall(kTRUE),
49  _nDim(0), _logfile(0),
50  _verbose(verbose)
51 {
52 
53  _evalCounter = 0 ;
54 
55  // Examine parameter list
56  RooArgSet* paramSet = _funct->getParameters(RooArgSet());
57  RooArgList paramList(*paramSet);
58  delete paramSet;
59 
60  _floatParamList = (RooArgList*) paramList.selectByAttrib("Constant",kFALSE);
61  if (_floatParamList->getSize()>1) {
62  _floatParamList->sort();
63  }
64  _floatParamList->setName("floatParamList");
65 
66  _constParamList = (RooArgList*) paramList.selectByAttrib("Constant",kTRUE);
67  if (_constParamList->getSize()>1) {
68  _constParamList->sort();
69  }
70  _constParamList->setName("constParamList");
71 
72  // Remove all non-RooRealVar parameters from list (MINUIT cannot handle them)
73  TIterator* pIter = _floatParamList->createIterator();
74  RooAbsArg* arg;
75  while ((arg=(RooAbsArg*)pIter->Next())) {
76  if (!arg->IsA()->InheritsFrom(RooAbsRealLValue::Class())) {
77  oocoutW(_context,Minimization) << "RooMinimizerFcn::RooMinimizerFcn: removing parameter "
78  << arg->GetName()
79  << " from list because it is not of type RooRealVar" << endl;
80  _floatParamList->remove(*arg);
81  }
82  }
83  delete pIter;
84 
85  _nDim = _floatParamList->getSize();
86 
87  updateFloatVec() ;
88 
89  // Save snapshot of initial lists
90  _initFloatParamList = (RooArgList*) _floatParamList->snapshot(kFALSE) ;
91  _initConstParamList = (RooArgList*) _constParamList->snapshot(kFALSE) ;
92 
93 }
94 
95 
96 
97 RooMinimizerFcn::RooMinimizerFcn(const RooMinimizerFcn& other) : ROOT::Math::IBaseFunctionMultiDim(other),
98  _evalCounter(other._evalCounter),
99  _funct(other._funct),
100  _context(other._context),
101  _maxFCN(other._maxFCN),
102  _numBadNLL(other._numBadNLL),
103  _printEvalErrors(other._printEvalErrors),
104  _doEvalErrorWall(other._doEvalErrorWall),
105  _nDim(other._nDim),
106  _logfile(other._logfile),
107  _verbose(other._verbose),
108  _floatParamVec(other._floatParamVec)
109 {
110  _floatParamList = new RooArgList(*other._floatParamList) ;
111  _constParamList = new RooArgList(*other._constParamList) ;
112  _initFloatParamList = (RooArgList*) other._initFloatParamList->snapshot(kFALSE) ;
113  _initConstParamList = (RooArgList*) other._initConstParamList->snapshot(kFALSE) ;
114 }
115 
116 
117 RooMinimizerFcn::~RooMinimizerFcn()
118 {
119  delete _floatParamList;
120  delete _initFloatParamList;
121  delete _constParamList;
122  delete _initConstParamList;
123 }
124 
125 
126 ROOT::Math::IBaseFunctionMultiDim* RooMinimizerFcn::Clone() const
127 {
128  return new RooMinimizerFcn(*this) ;
129 }
130 
131 
132 Bool_t RooMinimizerFcn::Synchronize(std::vector<ROOT::Fit::ParameterSettings>& parameters,
133  Bool_t optConst, Bool_t verbose)
134 {
135 
136  // Internal function to synchronize TMinimizer with current
137  // information in RooAbsReal function parameters
138 
139  Bool_t constValChange(kFALSE) ;
140  Bool_t constStatChange(kFALSE) ;
141 
142  Int_t index(0) ;
143 
144  // Handle eventual migrations from constParamList -> floatParamList
145  for(index= 0; index < _constParamList->getSize() ; index++) {
146 
147  RooRealVar *par= dynamic_cast<RooRealVar*>(_constParamList->at(index)) ;
148  if (!par) continue ;
149 
150  RooRealVar *oldpar= dynamic_cast<RooRealVar*>(_initConstParamList->at(index)) ;
151  if (!oldpar) continue ;
152 
153  // Test if constness changed
154  if (!par->isConstant()) {
155 
156  // Remove from constList, add to floatList
157  _constParamList->remove(*par) ;
158  _floatParamList->add(*par) ;
159  _initFloatParamList->addClone(*oldpar) ;
160  _initConstParamList->remove(*oldpar) ;
161  constStatChange=kTRUE ;
162  _nDim++ ;
163 
164  if (verbose) {
165  oocoutI(_context,Minimization) << "RooMinimizerFcn::synchronize: parameter "
166  << par->GetName() << " is now floating." << endl ;
167  }
168  }
169 
170  // Test if value changed
171  if (par->getVal()!= oldpar->getVal()) {
172  constValChange=kTRUE ;
173  if (verbose) {
174  oocoutI(_context,Minimization) << "RooMinimizerFcn::synchronize: value of constant parameter "
175  << par->GetName()
176  << " changed from " << oldpar->getVal() << " to "
177  << par->getVal() << endl ;
178  }
179  }
180 
181  }
182 
183  // Update reference list
184  *_initConstParamList = *_constParamList ;
185 
186  // Synchronize MINUIT with function state
187  // Handle floatParamList
188  for(index= 0; index < _floatParamList->getSize(); index++) {
189  RooRealVar *par= dynamic_cast<RooRealVar*>(_floatParamList->at(index)) ;
190 
191  if (!par) continue ;
192 
193  Double_t pstep(0) ;
194  Double_t pmin(0) ;
195  Double_t pmax(0) ;
196 
197  if(!par->isConstant()) {
198 
199  // Verify that floating parameter is indeed of type RooRealVar
200  if (!par->IsA()->InheritsFrom(RooRealVar::Class())) {
201  oocoutW(_context,Minimization) << "RooMinimizerFcn::fit: Error, non-constant parameter "
202  << par->GetName()
203  << " is not of type RooRealVar, skipping" << endl ;
204  _floatParamList->remove(*par);
205  index--;
206  _nDim--;
207  continue ;
208  }
209 
210  // Set the limits, if not infinite
211  if (par->hasMin() )
212  pmin = par->getMin();
213  if (par->hasMax() )
214  pmax = par->getMax();
215 
216  // Calculate step size
217  pstep = par->getError();
218  if(pstep <= 0) {
219  // Floating parameter without error estitimate
220  if (par->hasMin() && par->hasMax()) {
221  pstep= 0.1*(pmax-pmin);
222 
223  // Trim default choice of error if within 2 sigma of limit
224  if (pmax - par->getVal() < 2*pstep) {
225  pstep = (pmax - par->getVal())/2 ;
226  } else if (par->getVal() - pmin < 2*pstep) {
227  pstep = (par->getVal() - pmin )/2 ;
228  }
229 
230  // If trimming results in zero error, restore default
231  if (pstep==0) {
232  pstep= 0.1*(pmax-pmin);
233  }
234 
235  } else {
236  pstep=1 ;
237  }
238  if(verbose) {
239  oocoutW(_context,Minimization) << "RooMinimizerFcn::synchronize: WARNING: no initial error estimate available for "
240  << par->GetName() << ": using " << pstep << endl;
241  }
242  }
243  } else {
244  pmin = par->getVal() ;
245  pmax = par->getVal() ;
246  }
247 
248  // new parameter
249  if (index>=Int_t(parameters.size())) {
250 
251  if (par->hasMin() && par->hasMax()) {
252  parameters.push_back(ROOT::Fit::ParameterSettings(par->GetName(),
253  par->getVal(),
254  pstep,
255  pmin,pmax));
256  }
257  else {
258  parameters.push_back(ROOT::Fit::ParameterSettings(par->GetName(),
259  par->getVal(),
260  pstep));
261  if (par->hasMin() )
262  parameters.back().SetLowerLimit(pmin);
263  else if (par->hasMax() )
264  parameters.back().SetUpperLimit(pmax);
265  }
266 
267  continue;
268 
269  }
270 
271  Bool_t oldFixed = parameters[index].IsFixed();
272  Double_t oldVar = parameters[index].Value();
273  Double_t oldVerr = parameters[index].StepSize();
274  Double_t oldVlo = parameters[index].LowerLimit();
275  Double_t oldVhi = parameters[index].UpperLimit();
276 
277  if (par->isConstant() && !oldFixed) {
278 
279  // Parameter changes floating -> constant : update only value if necessary
280  if (oldVar!=par->getVal()) {
281  parameters[index].SetValue(par->getVal());
282  if (verbose) {
283  oocoutI(_context,Minimization) << "RooMinimizerFcn::synchronize: value of parameter "
284  << par->GetName() << " changed from " << oldVar
285  << " to " << par->getVal() << endl ;
286  }
287  }
288  parameters[index].Fix();
289  constStatChange=kTRUE ;
290  if (verbose) {
291  oocoutI(_context,Minimization) << "RooMinimizerFcn::synchronize: parameter "
292  << par->GetName() << " is now fixed." << endl ;
293  }
294 
295  } else if (par->isConstant() && oldFixed) {
296 
297  // Parameter changes constant -> constant : update only value if necessary
298  if (oldVar!=par->getVal()) {
299  parameters[index].SetValue(par->getVal());
300  constValChange=kTRUE ;
301 
302  if (verbose) {
303  oocoutI(_context,Minimization) << "RooMinimizerFcn::synchronize: value of fixed parameter "
304  << par->GetName() << " changed from " << oldVar
305  << " to " << par->getVal() << endl ;
306  }
307  }
308 
309  } else {
310  // Parameter changes constant -> floating
311  if (!par->isConstant() && oldFixed) {
312  parameters[index].Release();
313  constStatChange=kTRUE ;
314 
315  if (verbose) {
316  oocoutI(_context,Minimization) << "RooMinimizerFcn::synchronize: parameter "
317  << par->GetName() << " is now floating." << endl ;
318  }
319  }
320 
321  // Parameter changes constant -> floating : update all if necessary
322  if (oldVar!=par->getVal() || oldVlo!=pmin || oldVhi != pmax || oldVerr!=pstep) {
323  parameters[index].SetValue(par->getVal());
324  parameters[index].SetStepSize(pstep);
325  if (par->hasMin() && par->hasMax() )
326  parameters[index].SetLimits(pmin,pmax);
327  else if (par->hasMin() )
328  parameters[index].SetLowerLimit(pmin);
329  else if (par->hasMax() )
330  parameters[index].SetUpperLimit(pmax);
331  }
332 
333  // Inform user about changes in verbose mode
334  if (verbose) {
335  // if ierr<0, par was moved from the const list and a message was already printed
336 
337  if (oldVar!=par->getVal()) {
338  oocoutI(_context,Minimization) << "RooMinimizerFcn::synchronize: value of parameter "
339  << par->GetName() << " changed from " << oldVar << " to "
340  << par->getVal() << endl ;
341  }
342  if (oldVlo!=pmin || oldVhi!=pmax) {
343  oocoutI(_context,Minimization) << "RooMinimizerFcn::synchronize: limits of parameter "
344  << par->GetName() << " changed from [" << oldVlo << "," << oldVhi
345  << "] to [" << pmin << "," << pmax << "]" << endl ;
346  }
347 
348  // If oldVerr=0, then parameter was previously fixed
349  if (oldVerr!=pstep && oldVerr!=0) {
350  oocoutI(_context,Minimization) << "RooMinimizerFcn::synchronize: error/step size of parameter "
351  << par->GetName() << " changed from " << oldVerr << " to " << pstep << endl ;
352  }
353  }
354  }
355  }
356 
357  if (optConst) {
358  if (constStatChange) {
359 
360  RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
361 
362  oocoutI(_context,Minimization) << "RooMinimizerFcn::synchronize: set of constant parameters changed, rerunning const optimizer" << endl ;
363  _funct->constOptimizeTestStatistic(RooAbsArg::ConfigChange) ;
364  } else if (constValChange) {
365  oocoutI(_context,Minimization) << "RooMinimizerFcn::synchronize: constant parameter values changed, rerunning const optimizer" << endl ;
366  _funct->constOptimizeTestStatistic(RooAbsArg::ValueChange) ;
367  }
368 
369  RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
370 
371  }
372 
373  updateFloatVec() ;
374 
375  return 0 ;
376 
377 }
378 
379 Double_t RooMinimizerFcn::GetPdfParamVal(Int_t index)
380 {
381  // Access PDF parameter value by ordinal index (needed by MINUIT)
382 
383  return ((RooRealVar*)_floatParamList->at(index))->getVal() ;
384 }
385 
386 Double_t RooMinimizerFcn::GetPdfParamErr(Int_t index)
387 {
388  // Access PDF parameter error by ordinal index (needed by MINUIT)
389  return ((RooRealVar*)_floatParamList->at(index))->getError() ;
390 }
391 
392 
393 void RooMinimizerFcn::SetPdfParamErr(Int_t index, Double_t value)
394 {
395  // Modify PDF parameter error by ordinal index (needed by MINUIT)
396 
397  ((RooRealVar*)_floatParamList->at(index))->setError(value) ;
398 }
399 
400 
401 
402 void RooMinimizerFcn::ClearPdfParamAsymErr(Int_t index)
403 {
404  // Modify PDF parameter error by ordinal index (needed by MINUIT)
405 
406  ((RooRealVar*)_floatParamList->at(index))->removeAsymError() ;
407 }
408 
409 
410 void RooMinimizerFcn::SetPdfParamErr(Int_t index, Double_t loVal, Double_t hiVal)
411 {
412  // Modify PDF parameter error by ordinal index (needed by MINUIT)
413 
414  ((RooRealVar*)_floatParamList->at(index))->setAsymError(loVal,hiVal) ;
415 }
416 
417 
418 void RooMinimizerFcn::BackProp(const ROOT::Fit::FitResult &results)
419 {
420  // Transfer MINUIT fit results back into RooFit objects
421 
422  for (Int_t index= 0; index < _nDim; index++) {
423  Double_t value = results.Value(index);
424  SetPdfParamVal(index, value);
425 
426  // Set the parabolic error
427  Double_t err = results.Error(index);
428  SetPdfParamErr(index, err);
429 
430  Double_t eminus = results.LowerError(index);
431  Double_t eplus = results.UpperError(index);
432 
433  if(eplus > 0 || eminus < 0) {
434  // Store the asymmetric error, if it is available
435  SetPdfParamErr(index, eminus,eplus);
436  } else {
437  // Clear the asymmetric error
438  ClearPdfParamAsymErr(index) ;
439  }
440  }
441 
442 }
443 
444 Bool_t RooMinimizerFcn::SetLogFile(const char* inLogfile)
445 {
446  // Change the file name for logging of a RooMinimizer of all MINUIT steppings
447  // through the parameter space. If inLogfile is null, the current log file
448  // is closed and logging is stopped.
449 
450  if (_logfile) {
451  oocoutI(_context,Minimization) << "RooMinimizerFcn::setLogFile: closing previous log file" << endl ;
452  _logfile->close() ;
453  delete _logfile ;
454  _logfile = 0 ;
455  }
456  _logfile = new ofstream(inLogfile) ;
457  if (!_logfile->good()) {
458  oocoutI(_context,Minimization) << "RooMinimizerFcn::setLogFile: cannot open file " << inLogfile << endl ;
459  _logfile->close() ;
460  delete _logfile ;
461  _logfile= 0;
462  }
463 
464  return kFALSE ;
465 
466 }
467 
468 
469 void RooMinimizerFcn::ApplyCovarianceMatrix(TMatrixDSym& V)
470 {
471  // Apply results of given external covariance matrix. i.e. propagate its errors
472  // to all RRV parameter representations and give this matrix instead of the
473  // HESSE matrix at the next save() call
474 
475  for (Int_t i=0 ; i<_nDim ; i++) {
476  // Skip fixed parameters
477  if (_floatParamList->at(i)->isConstant()) {
478  continue ;
479  }
480  SetPdfParamErr(i, sqrt(V(i,i))) ;
481  }
482 
483 }
484 
485 
486 Bool_t RooMinimizerFcn::SetPdfParamVal(const Int_t &index, const Double_t &value) const
487 {
488  //RooRealVar* par = (RooRealVar*)_floatParamList->at(index);
489  RooRealVar* par = (RooRealVar*)_floatParamVec[index] ;
490 
491  if (par->getVal()!=value) {
492  if (_verbose) cout << par->GetName() << "=" << value << ", " ;
493 
494  par->setVal(value);
495  return kTRUE;
496  }
497 
498  return kFALSE;
499 }
500 
501 
502 
503 ////////////////////////////////////////////////////////////////////////////////
504 
505 void RooMinimizerFcn::updateFloatVec()
506 {
507  _floatParamVec.clear() ;
508  RooFIter iter = _floatParamList->fwdIterator() ;
509  RooAbsArg* arg ;
510  _floatParamVec = std::vector<RooAbsArg*>(_floatParamList->getSize()) ;
511  Int_t i(0) ;
512  while((arg=iter.next())) {
513  _floatParamVec[i++] = arg ;
514  }
515 }
516 
517 
518 
519 double RooMinimizerFcn::DoEval(const double *x) const
520 {
521 
522  // Set the parameter values for this iteration
523  for (int index = 0; index < _nDim; index++) {
524  if (_logfile) (*_logfile) << x[index] << " " ;
525  SetPdfParamVal(index,x[index]);
526  }
527 
528  // Calculate the function for these parameters
529  RooAbsReal::setHideOffset(kFALSE) ;
530  double fvalue = _funct->getVal();
531  RooAbsReal::setHideOffset(kTRUE) ;
532 
533  if (RooAbsReal::numEvalErrors()>0 || fvalue > 1e30) {
534 
535  if (_printEvalErrors>=0) {
536 
537  if (_doEvalErrorWall) {
538  oocoutW(_context,Minimization) << "RooMinimizerFcn: Minimized function has error status." << endl
539  << "Returning maximum FCN so far (" << _maxFCN
540  << ") to force MIGRAD to back out of this region. Error log follows" << endl ;
541  } else {
542  oocoutW(_context,Minimization) << "RooMinimizerFcn: Minimized function has error status but is ignored" << endl ;
543  }
544 
545  Bool_t first(kTRUE) ;
546  ooccoutW(_context,Minimization) << "Parameter values: " ;
547  for (const auto par : *_floatParamList) {
548  auto var = static_cast<const RooRealVar*>(par);
549  if (first) { first = kFALSE ; } else ooccoutW(_context,Minimization) << ", " ;
550  ooccoutW(_context,Minimization) << var->GetName() << "=" << var->getVal() ;
551  }
552  ooccoutW(_context,Minimization) << endl ;
553 
554  RooAbsReal::printEvalErrors(ooccoutW(_context,Minimization),_printEvalErrors) ;
555  ooccoutW(_context,Minimization) << endl ;
556  }
557 
558  if (_doEvalErrorWall) {
559  fvalue = _maxFCN+1;
560  }
561 
562  RooAbsReal::clearEvalErrorLog() ;
563  _numBadNLL++ ;
564  } else {
565  _maxFCN = std::max(fvalue, _maxFCN);
566  }
567 
568  // Optional logging
569  if (_logfile)
570  (*_logfile) << setprecision(15) << fvalue << setprecision(4) << endl;
571  if (_verbose) {
572  cout << "\nprevFCN" << (_funct->isOffsetting()?"-offset":"") << " = " << setprecision(10)
573  << fvalue << setprecision(4) << " " ;
574  cout.flush() ;
575  }
576 
577  _evalCounter++ ;
578 
579  return fvalue;
580 }
581 
582 #endif
583