RooAbsMinimizerFcn.cxx

Go to the documentation of this file.

/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 * @(#)root/roofitcore:$Id$
 * Authors:                                                                  *
 *   AL, Alfio Lazzaro,   INFN Milan,        alfio.lazzaro@mi.infn.it        *
 *   PB, Patrick Bos, Netherlands eScience Center, p.bos@esciencecenter.nl   *
 *                                                                           *
 *                                                                           *
 * Redistribution and use in source and binary forms,                        *
 * with or without modification, are permitted according to the terms        *
 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
 *****************************************************************************/
 
//////////////////////////////////////////////////////////////////////////////
//
// RooAbsMinimizerFcn is an interface class to the ROOT::Math function
// for minimization. It contains only the "logistics" of synchronizing
// between Minuit and RooFit. Its subclasses implement actual interfacing
// to Minuit by subclassing IMultiGenFunction or IMultiGradFunction.
//
 
#include "RooAbsMinimizerFcn.h"
 
#include "RooAbsArg.h"
#include "RooAbsPdf.h"
#include "RooArgSet.h"
#include "RooDataSet.h"
#include "RooRealVar.h"
#include "RooAbsRealLValue.h"
#include "RooMsgService.h"
#include "RooNaNPacker.h"
 
#include "TClass.h"
#include "TMatrixDSym.h"
 
#include <fstream>
#include <iomanip>
 
using std::endl, std::ofstream, std::cout;
 
RooAbsMinimizerFcn::RooAbsMinimizerFcn(RooArgList paramList, RooMinimizer *context) : _context(context)
{
   // Examine parameter list
   _floatParamList.reset(static_cast<RooArgList *>(paramList.selectByAttrib("Constant", false)));
   if (_floatParamList->size() > 1) {
      _floatParamList->sort();
   }
   _floatParamList->setName("floatParamList");
 
   _constParamList.reset(static_cast<RooArgList *>(paramList.selectByAttrib("Constant", true)));
   if (_constParamList->size() > 1) {
      _constParamList->sort();
   }
   _constParamList->setName("constParamList");
 
   // Remove all non-RooRealVar parameters from list (MINUIT cannot handle them)
   for (unsigned int i = 0; i < _floatParamList->size();) { // Note: Counting loop, since removing from collection!
      const RooAbsArg *arg = (*_floatParamList).at(i);
      if (!arg->IsA()->InheritsFrom(RooAbsRealLValue::Class())) {
         oocoutW(_context, Minimization) << "RooAbsMinimizerFcn::RooAbsMinimizerFcn: removing parameter "
                                         << arg->GetName() << " from list because it is not of type RooRealVar" << endl;
         _floatParamList->remove(*arg);
      } else {
         ++i;
      }
   }
 
   _nDim = _floatParamList->size();
 
   // Save snapshot of initial lists
   _initFloatParamList = std::make_unique<RooArgList>();
   _initConstParamList = std::make_unique<RooArgList>();
   _floatParamList->snapshot(*_initFloatParamList, false);
   _constParamList->snapshot(*_initConstParamList, false);
}
 
RooAbsMinimizerFcn::RooAbsMinimizerFcn(const RooAbsMinimizerFcn &other)
   : _context(other._context),
     _maxFCN(other._maxFCN),
     _funcOffset(other._funcOffset),
     _numBadNLL(other._numBadNLL),
     _evalCounter(other._evalCounter),
     _nDim(other._nDim),
     _optConst(other._optConst),
     _floatParamList(std::make_unique<RooArgList>(*other._floatParamList)),
     _constParamList(std::make_unique<RooArgList>(*other._constParamList)),
     _initFloatParamList(std::make_unique<RooArgList>()),
     _initConstParamList(std::make_unique<RooArgList>()),
     _logfile(other._logfile)
{
   other._initFloatParamList->snapshot(*_initFloatParamList, false);
   other._initConstParamList->snapshot(*_initConstParamList, false);
}
 
/// Internal function to synchronize TMinimizer with current
/// information in RooAbsReal function parameters
bool RooAbsMinimizerFcn::synchronizeParameterSettings(std::vector<ROOT::Fit::ParameterSettings> &parameters,
                                                      bool optConst)
{
   bool constValChange(false);
   bool constStatChange(false);
 
   // Handle eventual migrations from constParamList -> floatParamList
   for (std::size_t index = 0; index < _constParamList->size(); index++) {
 
      RooRealVar *par = dynamic_cast<RooRealVar *>(_constParamList->at(index));
      if (!par)
         continue;
 
      RooRealVar *oldpar = dynamic_cast<RooRealVar *>(_initConstParamList->at(index));
      if (!oldpar)
         continue;
 
      // Test if constness changed
      if (!par->isConstant()) {
 
         // Remove from constList, add to floatList
         _constParamList->remove(*par);
         _floatParamList->add(*par);
         _initFloatParamList->addClone(*oldpar);
         _initConstParamList->remove(*oldpar);
         constStatChange = true;
         _nDim++;
 
         if (cfg().verbose) {
            oocoutI(_context, Minimization)
               << "RooAbsMinimizerFcn::synchronize: parameter " << par->GetName() << " is now floating." << endl;
         }
      }
 
      // Test if value changed
      if (par->getVal() != oldpar->getVal()) {
         constValChange = true;
         if (cfg().verbose) {
            oocoutI(_context, Minimization)
               << "RooAbsMinimizerFcn::synchronize: value of constant parameter " << par->GetName() << " changed from "
               << oldpar->getVal() << " to " << par->getVal() << endl;
         }
      }
   }
 
   // Update reference list
   _initConstParamList->assign(*_constParamList);
 
   // Synchronize MINUIT with function state
   // Handle floatParamList
   for (std::size_t index = 0; index < _floatParamList->size(); index++) {
      RooRealVar *par = dynamic_cast<RooRealVar *>(_floatParamList->at(index));
 
      if (!par)
         continue;
 
      double pstep(0);
      double pmin(0);
      double pmax(0);
 
      if (!par->isConstant()) {
 
         // Verify that floating parameter is indeed of type RooRealVar
         if (!par->IsA()->InheritsFrom(RooRealVar::Class())) {
            oocoutW(_context, Minimization) << "RooAbsMinimizerFcn::fit: Error, non-constant parameter "
                                            << par->GetName() << " is not of type RooRealVar, skipping" << endl;
            _floatParamList->remove(*par);
            index--;
            _nDim--;
            continue;
         }
         // make sure the parameter are in dirty state to enable
         // a real NLL computation when the minimizer calls the function the first time
         // (see issue #7659)
         par->setValueDirty();
 
         // Set the limits, if not infinite
         if (par->hasMin())
            pmin = par->getMin();
         if (par->hasMax())
            pmax = par->getMax();
 
         // Calculate step size
         pstep = par->getError();
         if (pstep <= 0) {
            // Floating parameter without error estitimate
            if (par->hasMin() && par->hasMax()) {
               pstep = 0.1 * (pmax - pmin);
 
               // Trim default choice of error if within 2 sigma of limit
               if (pmax - par->getVal() < 2 * pstep) {
                  pstep = (pmax - par->getVal()) / 2;
               } else if (par->getVal() - pmin < 2 * pstep) {
                  pstep = (par->getVal() - pmin) / 2;
               }
 
               // If trimming results in zero error, restore default
               if (pstep == 0) {
                  pstep = 0.1 * (pmax - pmin);
               }
 
            } else {
               pstep = 1;
            }
            if (cfg().verbose) {
               oocoutW(_context, Minimization)
                  << "RooAbsMinimizerFcn::synchronize: WARNING: no initial error estimate available for "
                  << par->GetName() << ": using " << pstep << endl;
            }
         }
      } else {
         pmin = par->getVal();
         pmax = par->getVal();
      }
 
      // new parameter
      if (index >= parameters.size()) {
 
         if (par->hasMin() && par->hasMax()) {
            parameters.emplace_back(par->GetName(), par->getVal(), pstep, pmin, pmax);
         } else {
            parameters.emplace_back(par->GetName(), par->getVal(), pstep);
            if (par->hasMin()) {
               parameters.back().SetLowerLimit(pmin);
            } else if (par->hasMax()) {
               parameters.back().SetUpperLimit(pmax);
            }
         }
 
         continue;
      }
 
      bool oldFixed = parameters[index].IsFixed();
      double oldVar = parameters[index].Value();
      double oldVerr = parameters[index].StepSize();
      double oldVlo = parameters[index].LowerLimit();
      double oldVhi = parameters[index].UpperLimit();
 
      if (par->isConstant() && !oldFixed) {
 
         // Parameter changes floating -> constant : update only value if necessary
         if (oldVar != par->getVal()) {
            parameters[index].SetValue(par->getVal());
            if (cfg().verbose) {
               oocoutI(_context, Minimization)
                  << "RooAbsMinimizerFcn::synchronize: value of parameter " << par->GetName() << " changed from "
                  << oldVar << " to " << par->getVal() << endl;
            }
         }
         parameters[index].Fix();
         constStatChange = true;
         if (cfg().verbose) {
            oocoutI(_context, Minimization)
               << "RooAbsMinimizerFcn::synchronize: parameter " << par->GetName() << " is now fixed." << endl;
         }
 
      } else if (par->isConstant() && oldFixed) {
 
         // Parameter changes constant -> constant : update only value if necessary
         if (oldVar != par->getVal()) {
            parameters[index].SetValue(par->getVal());
            constValChange = true;
 
            if (cfg().verbose) {
               oocoutI(_context, Minimization)
                  << "RooAbsMinimizerFcn::synchronize: value of fixed parameter " << par->GetName() << " changed from "
                  << oldVar << " to " << par->getVal() << endl;
            }
         }
 
      } else {
         // Parameter changes constant -> floating
         if (!par->isConstant() && oldFixed) {
            parameters[index].Release();
            constStatChange = true;
 
            if (cfg().verbose) {
               oocoutI(_context, Minimization)
                  << "RooAbsMinimizerFcn::synchronize: parameter " << par->GetName() << " is now floating." << endl;
            }
         }
 
         // Parameter changes constant -> floating : update all if necessary
         if (oldVar != par->getVal() || oldVlo != pmin || oldVhi != pmax || oldVerr != pstep) {
            parameters[index].SetValue(par->getVal());
            parameters[index].SetStepSize(pstep);
            if (par->hasMin() && par->hasMax()) {
               parameters[index].SetLimits(pmin, pmax);
            } else if (par->hasMin()) {
               parameters[index].SetLowerLimit(pmin);
            } else if (par->hasMax()) {
               parameters[index].SetUpperLimit(pmax);
            }
         }
 
         // Inform user about changes in verbose mode
         if (cfg().verbose) {
            // if ierr<0, par was moved from the const list and a message was already printed
 
            if (oldVar != par->getVal()) {
               oocoutI(_context, Minimization)
                  << "RooAbsMinimizerFcn::synchronize: value of parameter " << par->GetName() << " changed from "
                  << oldVar << " to " << par->getVal() << endl;
            }
            if (oldVlo != pmin || oldVhi != pmax) {
               oocoutI(_context, Minimization)
                  << "RooAbsMinimizerFcn::synchronize: limits of parameter " << par->GetName() << " changed from ["
                  << oldVlo << "," << oldVhi << "] to [" << pmin << "," << pmax << "]" << endl;
            }
 
            // If oldVerr=0, then parameter was previously fixed
            if (oldVerr != pstep && oldVerr != 0) {
               oocoutI(_context, Minimization)
                  << "RooAbsMinimizerFcn::synchronize: error/step size of parameter " << par->GetName()
                  << " changed from " << oldVerr << " to " << pstep << endl;
            }
         }
      }
   }
 
   if (optConst) {
      optimizeConstantTerms(constStatChange, constValChange);
   }
 
   return false;
}
 
bool RooAbsMinimizerFcn::Synchronize(std::vector<ROOT::Fit::ParameterSettings> &parameters)
{
   return synchronizeParameterSettings(parameters, _optConst);
}
 
/// Modify PDF parameter error by ordinal index (needed by MINUIT)
void RooAbsMinimizerFcn::SetPdfParamErr(int index, double value)
{
   static_cast<RooRealVar *>(_floatParamList->at(index))->setError(value);
}
 
/// Modify PDF parameter error by ordinal index (needed by MINUIT)
void RooAbsMinimizerFcn::ClearPdfParamAsymErr(int index)
{
   static_cast<RooRealVar *>(_floatParamList->at(index))->removeAsymError();
}
 
/// Modify PDF parameter error by ordinal index (needed by MINUIT)
void RooAbsMinimizerFcn::SetPdfParamErr(int index, double loVal, double hiVal)
{
   static_cast<RooRealVar *>(_floatParamList->at(index))->setAsymError(loVal, hiVal);
}
 
/// Transfer MINUIT fit results back into RooFit objects.
void RooAbsMinimizerFcn::BackProp(const ROOT::Fit::FitResult &results)
{
   for (unsigned int index = 0; index < _nDim; index++) {
      double value = results.Value(index);
      SetPdfParamVal(index, value);
 
      // Set the parabolic error
      double err = results.Error(index);
      SetPdfParamErr(index, err);
 
      double eminus = results.LowerError(index);
      double eplus = results.UpperError(index);
 
      if (eplus > 0 || eminus < 0) {
         // Store the asymmetric error, if it is available
         SetPdfParamErr(index, eminus, eplus);
      } else {
         // Clear the asymmetric error
         ClearPdfParamAsymErr(index);
      }
   }
}
 
/// Change the file name for logging of a RooMinimizer of all MINUIT steppings
/// through the parameter space. If inLogfile is null, the current log file
/// is closed and logging is stopped.
bool RooAbsMinimizerFcn::SetLogFile(const char *inLogfile)
{
   if (_logfile) {
      oocoutI(_context, Minimization) << "RooAbsMinimizerFcn::setLogFile: closing previous log file" << endl;
      _logfile->close();
      delete _logfile;
      _logfile = nullptr;
   }
   _logfile = new ofstream(inLogfile);
   if (!_logfile->good()) {
      oocoutI(_context, Minimization) << "RooAbsMinimizerFcn::setLogFile: cannot open file " << inLogfile << endl;
      _logfile->close();
      delete _logfile;
      _logfile = nullptr;
   }
 
   return false;
}
 
/// Apply results of given external covariance matrix. i.e. propagate its errors
/// to all RRV parameter representations and give this matrix instead of the
/// HESSE matrix at the next save() call
void RooAbsMinimizerFcn::ApplyCovarianceMatrix(TMatrixDSym &V)
{
   for (unsigned int i = 0; i < _nDim; i++) {
      // Skip fixed parameters
      if (_floatParamList->at(i)->isConstant()) {
         continue;
      }
      SetPdfParamErr(i, sqrt(V(i, i)));
   }
}
 
/// Set value of parameter i.
bool RooAbsMinimizerFcn::SetPdfParamVal(int index, double value) const
{
   auto par = static_cast<RooRealVar *>(&(*_floatParamList)[index]);
 
   if (par->getVal() != value) {
      if (cfg().verbose)
         cout << par->GetName() << "=" << value << ", ";
 
      par->setVal(value);
      return true;
   }
 
   return false;
}
 
/// Print information about why evaluation failed.
/// Using _printEvalErrors, the number of errors printed can be steered.
/// Negative values disable printing.
void RooAbsMinimizerFcn::printEvalErrors() const
{
   if (cfg().printEvalErrors < 0)
      return;
 
   std::ostringstream msg;
   if (cfg().doEEWall) {
      msg << "RooAbsMinimizerFcn: Minimized function has error status." << endl
          << "Returning maximum FCN so far (" << _maxFCN
          << ") to force MIGRAD to back out of this region. Error log follows.\n";
   } else {
      msg << "RooAbsMinimizerFcn: Minimized function has error status but is ignored.\n";
   }
 
   msg << "Parameter values: ";
   for (const auto par : *_floatParamList) {
      auto var = static_cast<const RooRealVar *>(par);
      msg << "\t" << var->GetName() << "=" << var->getVal();
   }
   msg << std::endl;
 
   RooAbsReal::printEvalErrors(msg, cfg().printEvalErrors);
   ooccoutW(_context, Minimization) << msg.str() << endl;
}
 
/// Apply corrections on the fvalue if errors were signaled.
///
/// Two kinds of errors are possible: 1. infinite or nan values (the latter
/// can be a signaling nan, using RooNaNPacker) or 2. logEvalError-type errors.
/// Both are caught here and fvalue is updated so that Minuit in turn is nudged
/// to move the search outside of the problematic parameter space area.
double RooAbsMinimizerFcn::applyEvalErrorHandling(double fvalue) const
{
   if (!std::isfinite(fvalue) || RooAbsReal::numEvalErrors() > 0 || fvalue > 1e30) {
      printEvalErrors();
      RooAbsReal::clearEvalErrorLog();
      _numBadNLL++;
 
      if (cfg().doEEWall) {
         const double badness = RooNaNPacker::unpackNaN(fvalue);
         fvalue = (std::isfinite(_maxFCN) ? _maxFCN : 0.) + cfg().recoverFromNaN * badness;
      }
   } else {
      if (_evalCounter > 0 && _evalCounter == _numBadNLL) {
         // This is the first time we get a valid function value; while before, the
         // function was always invalid. For invalid  cases, we returned values > 0.
         // Now, we offset valid values such that they are < 0.
         _funcOffset = -fvalue;
      }
      fvalue += _funcOffset;
      _maxFCN = std::max(fvalue, _maxFCN);
   }
   return fvalue;
}
 
void RooAbsMinimizerFcn::finishDoEval() const
{
 
   if (_context->_loggingToDataSet) {
 
      if (!_context->_logDataSet) {
         const char *name = "minimizer_log_dataset";
         _context->_logDataSet = std::make_unique<RooDataSet>(name, name, *_floatParamList);
      }
 
      _context->_logDataSet->add(*_floatParamList);
   }
 
   _evalCounter++;
}
 
void RooAbsMinimizerFcn::setOptimizeConst(int flag)
{
   auto ctx = _context->makeEvalErrorContext();
 
   if (_optConst && !flag) {
      if (_context->getPrintLevel() > -1) {
         oocoutI(_context, Minimization) << "RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization"
                                         << endl;
      }
      setOptimizeConstOnFunction(RooAbsArg::DeActivate, true);
      _optConst = flag;
   } else if (!_optConst && flag) {
      if (_context->getPrintLevel() > -1) {
         oocoutI(_context, Minimization) << "RooAbsMinimizerFcn::setOptimizeConst: activating const optimization"
                                         << endl;
      }
      setOptimizeConstOnFunction(RooAbsArg::Activate, flag > 1);
      _optConst = flag;
   } else if (_optConst && flag) {
      if (_context->getPrintLevel() > -1) {
         oocoutI(_context, Minimization) << "RooAbsMinimizerFcn::setOptimizeConst: const optimization already active"
                                         << endl;
      }
   } else {
      if (_context->getPrintLevel() > -1) {
         oocoutI(_context, Minimization) << "RooAbsMinimizerFcn::setOptimizeConst: const optimization wasn't active"
                                         << endl;
      }
   }
}
 
void RooAbsMinimizerFcn::optimizeConstantTerms(bool constStatChange, bool constValChange)
{
   auto ctx = _context->makeEvalErrorContext();
 
   if (constStatChange) {
 
      oocoutI(_context, Minimization)
         << "RooAbsMinimizerFcn::optimizeConstantTerms: set of constant parameters changed, rerunning const optimizer"
         << endl;
      setOptimizeConstOnFunction(RooAbsArg::ConfigChange, true);
   } else if (constValChange) {
      oocoutI(_context, Minimization)
         << "RooAbsMinimizerFcn::optimizeConstantTerms: constant parameter values changed, rerunning const optimizer"
         << endl;
      setOptimizeConstOnFunction(RooAbsArg::ValueChange, true);
   }
}
 
std::vector<double> RooAbsMinimizerFcn::getParameterValues() const
{
   // TODO: make a cache for this somewhere so it doesn't have to be recreated on each call
   std::vector<double> values;
   values.reserve(_nDim);
 
   for (std::size_t index = 0; index < _nDim; ++index) {
      RooRealVar *par = static_cast<RooRealVar *>(_floatParamList->at(index));
      values.push_back(par->getVal());
   }
 
   return values;
}