From dd51352075abc5b38ee13191713751e56696f23a Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Thu, 5 Oct 2023 13:55:26 +0200 Subject: [PATCH] [RF][HF] Clearly mark interpolation code 3 as unknown As noted in GitHub issue #7103, the interpolation code 3 is the same as code 2 in the `FlexibleInterpVar` and the `PiecewiseInterpolation` classes. According to some comments in the source code, interpolation code 3 was meant to be "a parabolic version of log-normal". There were two options to fix this: 1) Actually implement this parabolic interpolation with linear extrapolation, analogous to code 2 but in log space. 2) Clearly mark interpolation code 3 as non-existing. This commit implements solution 2, because the interpolation code 3 is not mentioned anywhere outside the ROOT source code. Especially not is the HistFactory paper: https://cds.cern.ch/record/1456844/files/CERN-OPEN-2012-016.pdf Implementing a new interpolation scheme that apparently was not needed in the last 10 years anyway would increase the burden on the user to understand all these different settings unnecessarily. Closes #7103. --- .../RooStats/HistFactory/FlexibleInterpVar.h | 4 + .../HistFactory/PiecewiseInterpolation.h | 6 +- roofit/histfactory/src/FlexibleInterpVar.cxx | 77 +++++++++-------- .../src/HistoToWorkspaceFactoryFast.cxx | 2 +- .../src/PiecewiseInterpolation.cxx | 86 +++++++++---------- .../roofitcore/inc/RooFit/Detail/MathFuncs.h | 16 +--- 6 files changed, 95 insertions(+), 96 deletions(-) diff --git a/roofit/histfactory/inc/RooStats/HistFactory/FlexibleInterpVar.h b/roofit/histfactory/inc/RooStats/HistFactory/FlexibleInterpVar.h index 7fbdf1ddacb30..46285e2f60b68 100644 --- a/roofit/histfactory/inc/RooStats/HistFactory/FlexibleInterpVar.h +++ b/roofit/histfactory/inc/RooStats/HistFactory/FlexibleInterpVar.h @@ -71,6 +71,10 @@ namespace HistFactory{ double evaluate() const override; + private: + + void setInterpCodeForParam(int iParam, int code); + ClassDefOverride(RooStats::HistFactory::FlexibleInterpVar,2); // flexible interpolation }; } diff --git a/roofit/histfactory/inc/RooStats/HistFactory/PiecewiseInterpolation.h b/roofit/histfactory/inc/RooStats/HistFactory/PiecewiseInterpolation.h index cd47e4dcd12c4..5676b7f4d1c1d 100644 --- a/roofit/histfactory/inc/RooStats/HistFactory/PiecewiseInterpolation.h +++ b/roofit/histfactory/inc/RooStats/HistFactory/PiecewiseInterpolation.h @@ -61,7 +61,7 @@ class PiecewiseInterpolation : public RooAbsReal { void setPositiveDefinite(bool flag=true){_positiveDefinite=flag;} bool positiveDefinite() const {return _positiveDefinite;} - void setInterpCode(RooAbsReal& param, int code, bool silent=false); + void setInterpCode(RooAbsReal& param, int code, bool silent=true); void setAllInterpCodes(int code); void printAllInterpCodes(); @@ -102,6 +102,10 @@ class PiecewiseInterpolation : public RooAbsReal { double evaluate() const override; void doEval(RooFit::EvalContext &) const override; +private: + + void setInterpCodeForParam(int iParam, int code); + ClassDefOverride(PiecewiseInterpolation,4) // Sum of RooAbsReal objects }; diff --git a/roofit/histfactory/src/FlexibleInterpVar.cxx b/roofit/histfactory/src/FlexibleInterpVar.cxx index 11c97bdacec02..d65dfd2e64d01 100644 --- a/roofit/histfactory/src/FlexibleInterpVar.cxx +++ b/roofit/histfactory/src/FlexibleInterpVar.cxx @@ -54,10 +54,10 @@ FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title, FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title, const RooArgList& paramList, double argNominal, std::vector const& lowVec, std::vector const& highVec, - std::vector const& code) : + std::vector const& codes) : RooAbsReal(name, title), _paramList("paramList","List of paramficients",this), - _nominal(argNominal), _low(lowVec), _high(highVec), _interpCode(code) + _nominal(argNominal), _low(lowVec), _high(highVec) { for (auto param : paramList) { if (!dynamic_cast(param)) { @@ -69,6 +69,11 @@ FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title, _paramList.add(*param) ; } + _interpCode.resize(_paramList.size()); + for (std::size_t i = 0; i < codes.size(); ++i) { + setInterpCodeForParam(i, codes[i]); + } + if (_low.size() != _paramList.size() || _low.size() != _high.size() || _low.size() != _interpCode.size()) { coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") invalid input std::vectors " << std::endl; R__ASSERT(_low.size() == _paramList.size()); @@ -109,31 +114,43 @@ FlexibleInterpVar::~FlexibleInterpVar() TRACE_DESTROY; } - -//////////////////////////////////////////////////////////////////////////////// - -void FlexibleInterpVar::setInterpCode(RooAbsReal& param, int code){ - int index = _paramList.index(¶m); - if(index<0){ - coutE(InputArguments) << "FlexibleInterpVar::setInterpCode ERROR: " << param.GetName() - << " is not in list" << std::endl; - } else if(_interpCode.at(index) != code){ - coutI(InputArguments) << "FlexibleInterpVar::setInterpCode : " << param.GetName() - << " is now " << code << std::endl; - _interpCode.at(index) = code; - // GHL: Adding suggestion by Swagato: - setValueDirty(); - } +void FlexibleInterpVar::setInterpCode(RooAbsReal ¶m, int code) +{ + int index = _paramList.index(¶m); + if (index < 0) { + coutE(InputArguments) << "FlexibleInterpVar::setInterpCode ERROR: " << param.GetName() << " is not in list" + << std::endl; + return; + } + setInterpCodeForParam(index, code); } -//////////////////////////////////////////////////////////////////////////////// +void FlexibleInterpVar::setAllInterpCodes(int code) +{ + for (std::size_t i = 0; i < _interpCode.size(); ++i) { + setInterpCodeForParam(i, code); + } +} -void FlexibleInterpVar::setAllInterpCodes(int code){ - for(unsigned int i=0; i<_interpCode.size(); ++i){ - _interpCode.at(i) = code; - } - // GHL: Adding suggestion by Swagato: - setValueDirty(); +void FlexibleInterpVar::setInterpCodeForParam(int iParam, int code) +{ + RooAbsArg const ¶m = _paramList[iParam]; + if (code < 0 || code > 5) { + coutE(InputArguments) << "FlexibleInterpVar::setInterpCode ERROR: " << param.GetName() + << " with unknown interpolation code " << code << ", keeping current code " + << _interpCode[iParam] << std::endl; + return; + } + if (code == 3) { + // In the past, code 3 was equivalent to code 2, which confused users. + // Now, we just say that code 3 doesn't exist and default to code 2 in + // that case for backwards compatible behavior. + coutE(InputArguments) << "FlexibleInterpVar::setInterpCode ERROR: " << param.GetName() + << " with unknown interpolation code " << code << ", defaulting to code 2" << std::endl; + code = 2; + } + _interpCode.at(iParam) = code; + setValueDirty(); } //////////////////////////////////////////////////////////////////////////////// @@ -198,10 +215,6 @@ double FlexibleInterpVar::evaluate() const double total(_nominal); for (std::size_t i = 0; i < _paramList.size(); ++i) { int code = _interpCode[i]; - if (code < 0 || code > 4) { - coutE(InputArguments) << "FlexibleInterpVar::evaluate ERROR: param " << i - << " with unknown interpolation code" << std::endl; - } // To get consistent codes with the PiecewiseInterpolation if (code == 4) { code = 5; @@ -223,10 +236,6 @@ void FlexibleInterpVar::translate(RooFit::Detail::CodeSquashContext &ctx) const unsigned int n = _interpCode.size(); int interpCode = _interpCode[0]; - if (interpCode < 0 || interpCode > 4) { - coutE(InputArguments) << "FlexibleInterpVar::evaluate ERROR: param " << 0 - << " with unknown interpolation code" << std::endl; - } // To get consistent codes with the PiecewiseInterpolation if (interpCode == 4) { interpCode = 5; @@ -251,10 +260,6 @@ void FlexibleInterpVar::doEval(RooFit::EvalContext &ctx) const for (std::size_t i = 0; i < _paramList.size(); ++i) { int code = _interpCode[i]; - if (code < 0 || code > 4) { - coutE(InputArguments) << "FlexibleInterpVar::evaluate ERROR: param " << i - << " with unknown interpolation code" << std::endl; - } // To get consistent codes with the PiecewiseInterpolation if (code == 4) { code = 5; diff --git a/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx b/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx index ef2eb70aacd18..d14eda1391097 100644 --- a/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx +++ b/roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx @@ -603,7 +603,7 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo assert(lowVec.size() == params.size()); FlexibleInterpVar interp( (interpName).c_str(), "", params, 1., lowVec, highVec); - interp.setAllInterpCodes(4); // LM: change to 4 (piece-wise linear to 6th order polynomial interpolation + linear extrapolation ) + interp.setAllInterpCodes(4); // LM: change to 4 (piece-wise exponential to 6th order polynomial interpolation + exponential extrapolation ) //interp.setAllInterpCodes(0); // simple linear interpolation proto.import(interp); // params have already been imported in first loop of this function } else{ diff --git a/roofit/histfactory/src/PiecewiseInterpolation.cxx b/roofit/histfactory/src/PiecewiseInterpolation.cxx index 2421dc2767674..3dca4123f8e4f 100644 --- a/roofit/histfactory/src/PiecewiseInterpolation.cxx +++ b/roofit/histfactory/src/PiecewiseInterpolation.cxx @@ -161,14 +161,8 @@ double PiecewiseInterpolation::evaluate() const auto param = static_cast(_paramSet.at(i)); auto low = static_cast(_lowSet.at(i)); auto high = static_cast(_highSet.at(i)); - Int_t icode = _interpCode[i] ; - - if(icode < 0 || icode > 5) { - coutE(InputArguments) << "PiecewiseInterpolation::evaluate ERROR: " << param->GetName() - << " with unknown interpolation code" << icode << endl ; - } using RooFit::Detail::MathFuncs::flexibleInterpSingle; - sum += flexibleInterpSingle(icode, low->getVal(), high->getVal(), 1.0, nominal, param->getVal(), sum); + sum += flexibleInterpSingle(_interpCode[i], low->getVal(), high->getVal(), 1.0, nominal, param->getVal(), sum); } if(_positiveDefinite && (sum<0)){ @@ -190,10 +184,6 @@ void PiecewiseInterpolation::translate(RooFit::Detail::CodeSquashContext &ctx) c std::string resName = "total_" + ctx.getTmpVarName(); for (std::size_t i = 0; i < n; ++i) { - if (_interpCode[i] < 0 || _interpCode[i] > 5) { - coutE(InputArguments) << "PiecewiseInterpolation::evaluate ERROR: " << _paramSet[i].GetName() - << " with unknown interpolation code" << _interpCode[i] << endl; - } if (_interpCode[i] != _interpCode[0]) { coutE(InputArguments) << "FlexibleInterpVar::evaluate ERROR: Code Squashing AD does not yet support having " "different interpolation codes for the same class object " @@ -277,18 +267,10 @@ void PiecewiseInterpolation::doEval(RooFit::EvalContext &ctx) const auto param = ctx.at(_paramSet.at(i)); auto low = ctx.at(_lowSet.at(i)); auto high = ctx.at(_highSet.at(i)); - const int icode = _interpCode[i]; - - if (icode < 0 || icode > 5) { - coutE(InputArguments) << "PiecewiseInterpolation::doEval(): " << _paramSet[i].GetName() - << " with unknown interpolation code" << icode << std::endl; - throw std::invalid_argument("PiecewiseInterpolation::doEval() got invalid interpolation code " + - std::to_string(icode)); - } for (std::size_t j = 0; j < sum.size(); ++j) { using RooFit::Detail::MathFuncs::flexibleInterpSingle; - sum[j] += flexibleInterpSingle(icode, broadcast(low, j), broadcast(high, j), 1.0, broadcast(nominal, j), + sum[j] += flexibleInterpSingle(_interpCode[i], broadcast(low, j), broadcast(high, j), 1.0, broadcast(nominal, j), broadcast(param, j), sum[j]); } } @@ -347,7 +329,7 @@ Int_t PiecewiseInterpolation::getAnalyticalIntegralWN(RooArgSet& allVars, RooArg // KC: check if interCode=0 for all - for (auto it = _paramSet.begin(); it != _paramSet.end(); ++it) { + for (auto it = _paramSet.begin(); it != _paramSet.end(); ++it) { if (!_interpCode.empty() && _interpCode[it - _paramSet.begin()] != 0) { // can't factorize integral cout << "can't factorize integral" << endl; @@ -371,7 +353,7 @@ Int_t PiecewiseInterpolation::getAnalyticalIntegralWN(RooArgSet& allVars, RooArg // Make list of function projection and normalization integrals RooAbsReal *func ; - // do variations + // do variations for (auto it = _paramSet.begin(); it != _paramSet.end(); ++it) { auto i = it - _paramSet.begin(); @@ -491,12 +473,12 @@ double PiecewiseInterpolation::analyticalIntegralWN(Int_t code, const RooArgSet* // now get low/high variations // KC: old interp code with new iterator - + i = 0; for (auto const *param : static_range_cast(_paramSet)) { low = static_cast(cache->_lowIntList.at(i)); high = static_cast(cache->_highIntList.at(i)); - + if(param->getVal() > 0) { value += param->getVal()*(high->getVal() - nominal); } else { @@ -573,32 +555,44 @@ double PiecewiseInterpolation::analyticalIntegralWN(Int_t code, const RooArgSet* return value; } - -//////////////////////////////////////////////////////////////////////////////// - -void PiecewiseInterpolation::setInterpCode(RooAbsReal& param, int code, bool silent){ - int index = _paramSet.index(¶m); - if(index<0){ - coutE(InputArguments) << "PiecewiseInterpolation::setInterpCode ERROR: " << param.GetName() - << " is not in list" << endl ; - } else { - if(!silent){ - coutW(InputArguments) << "PiecewiseInterpolation::setInterpCode : " << param.GetName() - << " is now " << code << endl ; - } - _interpCode.at(index) = code; - } +void PiecewiseInterpolation::setInterpCode(RooAbsReal ¶m, int code, bool /*silent*/) +{ + int index = _paramSet.index(¶m); + if (index < 0) { + coutE(InputArguments) << "PiecewiseInterpolation::setInterpCode ERROR: " << param.GetName() << " is not in list" + << std::endl; + return; + } + setInterpCodeForParam(index, code); } - -//////////////////////////////////////////////////////////////////////////////// - -void PiecewiseInterpolation::setAllInterpCodes(int code){ - for(unsigned int i=0; i<_interpCode.size(); ++i){ - _interpCode.at(i) = code; - } +void PiecewiseInterpolation::setAllInterpCodes(int code) +{ + for (std::size_t i = 0; i < _interpCode.size(); ++i) { + setInterpCodeForParam(i, code); + } } +void PiecewiseInterpolation::setInterpCodeForParam(int iParam, int code) +{ + RooAbsArg const ¶m = _paramSet[iParam]; + if (code < 0 || code > 5) { + coutE(InputArguments) << "PiecewiseInterpolation::setInterpCode ERROR: " << param.GetName() + << " with unknown interpolation code " << code << ", keeping current code " + << _interpCode[iParam] << std::endl; + return; + } + if (code == 3) { + // In the past, code 3 was equivalent to code 2, which confused users. + // Now, we just say that code 3 doesn't exist and default to code 2 in + // that case for backwards compatible behavior. + coutE(InputArguments) << "PiecewiseInterpolation::setInterpCode ERROR: " << param.GetName() + << " with unknown interpolation code " << code << ", defaulting to code 2" << std::endl; + code = 2; + } + _interpCode.at(iParam) = code; + setValueDirty(); +} //////////////////////////////////////////////////////////////////////////////// diff --git a/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h b/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h index b68cfa958c781..625349b37ef08 100644 --- a/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h +++ b/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h @@ -239,18 +239,10 @@ inline double flexibleInterpSingle(unsigned int code, double low, double high, d } else { return a * std::pow(paramVal, 2) + b * paramVal + c; } - } else if (code == 3) { - // parabolic version of log-normal - double a = 0.5 * (high + low) - nominal; - double b = 0.5 * (high - low); - double c = 0; - if (paramVal > 1) { - return (2 * a + b) * (paramVal - 1) + high - nominal; - } else if (paramVal < -1) { - return -1 * (2 * a - b) * (paramVal + 1) + low - nominal; - } else { - return a * std::pow(paramVal, 2) + b * paramVal + c; - } + // According to an old comment in the source code, code 3 was apparently + // meant to be a "parabolic version of log-normal", but it never got + // implemented. If someone would need it, it could be implemented as doing + // code 2 in log space. } else if (code == 4) { double x = paramVal; if (x >= boundary) {