diff --git a/interface/CMSHistErrorPropagator.h b/interface/CMSHistErrorPropagator.h index e709adf0b56..08835c5c5c0 100644 --- a/interface/CMSHistErrorPropagator.h +++ b/interface/CMSHistErrorPropagator.h @@ -71,7 +71,8 @@ class CMSHistErrorPropagator : public RooAbsReal { RooArgList wrapperList() const; RooArgList const& coefList() const { return coeffs_; } RooArgList const& funcList() const { return funcs_; } - + RooAbsReal const& getXVar() const { return *x_; } + std::map getProcessNorms() const; friend class CMSHistV; diff --git a/interface/CMSHistSum.h b/interface/CMSHistSum.h index 082b9264183..f7cf707499e 100644 --- a/interface/CMSHistSum.h +++ b/interface/CMSHistSum.h @@ -74,6 +74,8 @@ class CMSHistSum : public RooAbsReal { RooAbsReal const& getXVar() const { return x_.arg(); } + std::vector getFuncValList(std::size_t fnIdx); + static void EnableFastVertical(); friend class CMSHistV; diff --git a/interface/CombineCodegenImpl.h b/interface/CombineCodegenImpl.h index 49ac1882337..5cdeebe2aaf 100644 --- a/interface/CombineCodegenImpl.h +++ b/interface/CombineCodegenImpl.h @@ -4,6 +4,9 @@ #include class AsymPow; +class CMSHistErrorPropagator; +class CMSHistFunc; +class CMSHistSum; class FastVerticalInterpHistPdf2; class FastVerticalInterpHistPdf2D2; class ProcessNormalization; @@ -14,11 +17,17 @@ namespace RooFit::Experimental { class CodegenContext; void codegenImpl(AsymPow& arg, CodegenContext& ctx); + void codegenImpl(CMSHistErrorPropagator& arg, CodegenContext& ctx); + void codegenImpl(CMSHistFunc& arg, CodegenContext& ctx); + void codegenImpl(CMSHistSum& arg, CodegenContext& ctx); void codegenImpl(FastVerticalInterpHistPdf2& arg, CodegenContext& ctx); void codegenImpl(FastVerticalInterpHistPdf2D2& arg, CodegenContext& ctx); void codegenImpl(ProcessNormalization& arg, CodegenContext& ctx); void codegenImpl(VerticalInterpPdf& arg, CodegenContext& ctx); + std::string codegenIntegralImpl(CMSHistErrorPropagator& arg, int code, const char* rangeName, CodegenContext& ctx); + std::string codegenIntegralImpl(CMSHistFunc& arg, int code, const char* rangeName, CodegenContext& ctx); + std::string codegenIntegralImpl(CMSHistSum& arg, int code, const char* rangeName, CodegenContext& ctx); std::string codegenIntegralImpl(VerticalInterpPdf& arg, int code, const char* rangeName, CodegenContext& ctx); } // namespace RooFit::Experimental diff --git a/interface/CombineMathFuncs.h b/interface/CombineMathFuncs.h index 0dddea890e4..b7508dbbdc3 100644 --- a/interface/CombineMathFuncs.h +++ b/interface/CombineMathFuncs.h @@ -7,6 +7,8 @@ #include #include +#include + #include namespace RooFit { @@ -306,6 +308,32 @@ inline Double_t verticalInterpPdfIntegral(double const* coefList, std::size_t nC return result > 0. ? result : integralFloorVal; } +inline double cmsHistFunc(double x, std::size_t nBins, double const* binEdges, double const *values) { + // I guess a "CMS hist func" is just looking up some values in some bin? + unsigned int binIdx = RooFit::Detail::MathFuncs::binNumber(x, 1.0, binEdges, nBins + 1, nBins, 0); + return values[binIdx]; +} + +inline double cmsHistErrorPropagator(double x, std::size_t nFuncs, double const* coefList, double const* funcList) { + // My naive understanding of the logic: multiply functions with coefficients and sum up + double out = 0.; + for (std::size_t i = 0; i < nFuncs; ++i) { + out += coefList[i] * funcList[i]; + } + return out; +} + +inline double cmsHistSum( + double x, std::size_t nBins, std::size_t nSamples, double* coeffs, double const* binEdges, double const* values) { + unsigned int binIdx = RooFit::Detail::MathFuncs::binNumber(x, 1.0, binEdges, nBins + 1, nBins, 0); + double val = 0.; + for (std::size_t iSample = 0; iSample < nSamples; ++iSample) + val += coeffs[iSample] * values[iSample * nBins + binIdx]; + + return val; +} + + } // namespace MathFuncs } // namespace Detail } // namespace RooFit diff --git a/interface/FastTemplate_Old.h b/interface/FastTemplate_Old.h index 5461d8f2109..6f5b86f51bc 100644 --- a/interface/FastTemplate_Old.h +++ b/interface/FastTemplate_Old.h @@ -76,13 +76,15 @@ class FastTemplate { void Dump() const ; + AT const &GetValues() const { return values_; } + protected: unsigned int size_; AT values_; }; class FastHisto : public FastTemplate { public: - FastHisto() : FastTemplate(), binEdges_(), binWidths_() {} + FastHisto() = default; FastHisto(const TH1 &hist) ; FastHisto(const FastHisto &other) ; FastHisto & operator=(const FastHisto &other) { @@ -124,6 +126,8 @@ class FastHisto : public FastTemplate { const T & GetEdge(unsigned int i) const { return binEdges_[i]; } const T & GetWidth(unsigned int i) const { return binWidths_[i]; } + AT const &GetBinEdges() const { return binEdges_; } + private: AT binEdges_; AT binWidths_; diff --git a/src/CMSHistSum.cc b/src/CMSHistSum.cc index 9d46a78fa1f..41f9001ef6c 100644 --- a/src/CMSHistSum.cc +++ b/src/CMSHistSum.cc @@ -412,6 +412,24 @@ void CMSHistSum::updateCache() const { } } +std::vector CMSHistSum::getFuncValList(std::size_t fnIdx) { + staging_ = compcache_[fnIdx]; + if (vtype_[fnIdx] == CMSHistFunc::VerticalSetting::LogQuadLinear) { + staging_.Exp(); + staging_.Scale(storage_[process_fields_[fnIdx]].Integral() / staging_.Integral()); + } + staging_.CropUnderflows(); + std::vector result = staging_.GetValues(); + for (unsigned j = 0; j < bintypes_.size(); ++j) { + if (bintypes_[j][0] == 1) { + double x = vbinpars_[j][0]->getVal(); + result[j] += binerrors_[fnIdx][j] * x; + } else if (bintypes_[j][0] == 2 || bintypes_[j][0] == 3) + result[j] += scaledbinmods_[fnIdx][j]; + } + return result; +} + void CMSHistSum::runBarlowBeeston() const { updateCache(); const unsigned n = bb_.use.size(); diff --git a/src/CombineCodegenImpl.cxx b/src/CombineCodegenImpl.cxx index 0d1e0c09f3f..922fb3e0c77 100644 --- a/src/CombineCodegenImpl.cxx +++ b/src/CombineCodegenImpl.cxx @@ -5,10 +5,13 @@ #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 36, 0) #include "../interface/AsymPow.h" +#include "../interface/CMSHistErrorPropagator.h" +#include "../interface/CMSHistFunc.h" +#include "../interface/CMSHistSum.h" +#include "../interface/CombineMathFuncs.h" #include "../interface/ProcessNormalization.h" #include "../interface/VerticalInterpHistPdf.h" #include "../interface/VerticalInterpPdf.h" -#include "../interface/CombineMathFuncs.h" #include @@ -257,4 +260,68 @@ std::string RooFit::Experimental::codegenIntegralImpl(VerticalInterpPdf& arg, arg.quadraticAlgo()); } +void RooFit::Experimental::codegenImpl(CMSHistFunc& arg, CodegenContext& ctx) { + arg.evaluate(); // trigger cache() creation + std::vector const& edges = arg.cache().GetBinEdges(); + + // I don't know if these values are actually constant and we can take them + // here to hardcode into the generated code... + auto const& values = arg.cache().GetValues(); + + ctx.addResult(&arg, + ctx.buildCall("RooFit::Detail::MathFuncs::cmsHistFunc", + arg.getXVar(), + edges.size() - 1, + edges, + values + )); +} + +void RooFit::Experimental::codegenImpl(CMSHistErrorPropagator& arg, CodegenContext& ctx) { + ctx.addResult(&arg, + ctx.buildCall("RooFit::Detail::MathFuncs::cmsHistErrorPropagator", + arg.getXVar(), + arg.coefList().size(), + arg.coefList(), + arg.funcList())); +} + +std::string RooFit::Experimental::codegenIntegralImpl(CMSHistErrorPropagator& arg, + int code, + const char* rangeName, + CodegenContext& ctx) { + return "2.0"; // TODO: dummy for now +} + +void RooFit::Experimental::codegenImpl(CMSHistSum& arg, CodegenContext& ctx) { + arg.evaluate(); + std::vector const& edges = arg.cache().GetBinEdges(); + std::size_t nBins = edges.size() - 1; + RooArgList const& coefs = arg.coefList(); + std::size_t nSamples = coefs.size(); + std::vector values(nBins * nSamples); + for (std::size_t iSamples = 0; iSamples < nSamples; ++iSamples) { + std::vector sampleValues = arg.getFuncValList(iSamples); + for (std::size_t iBin = 0; iBin < nBins; ++iBin) + values[iBin + iSamples * nBins] = sampleValues[iBin]; + } + + ctx.addResult(&arg, + ctx.buildCall("RooFit::Detail::MathFuncs::cmsHistSum", + arg.getXVar(), + nBins, + nSamples, + coefs, + edges, + values + )); +} + +std::string RooFit::Experimental::codegenIntegralImpl(CMSHistSum& arg, + int code, + const char* rangeName, + CodegenContext& ctx) { + return "3.0"; // TODO: dummy for now +} + #endif // ROOT_VERSION_CODE >= ROOT_VERSION(6,36,0) diff --git a/src/FastTemplate_Old.cc b/src/FastTemplate_Old.cc index 7d617ea9b60..682d4e09165 100644 --- a/src/FastTemplate_Old.cc +++ b/src/FastTemplate_Old.cc @@ -73,7 +73,7 @@ FastHisto::FastHisto(const FastHisto &other) : } int FastHisto::FindBin(const T &x) const { - auto match = std::lower_bound(binEdges_.begin(), binEdges_.end(), x); + auto match = std::upper_bound(binEdges_.begin(), binEdges_.end(), x); if (match == binEdges_.begin()) return -1; if (match == binEdges_.end()) return values_.size(); return match - binEdges_.begin() - 1; @@ -81,7 +81,7 @@ int FastHisto::FindBin(const T &x) const { FastHisto::T FastHisto::GetAt(const T &x) const { - auto match = std::lower_bound(binEdges_.begin(), binEdges_.end(), x); + auto match = std::upper_bound(binEdges_.begin(), binEdges_.end(), x); if (match == binEdges_.begin() || match == binEdges_.end()) return T(0.0); return values_[match - binEdges_.begin() - 1]; } @@ -140,10 +140,10 @@ FastHisto2D::FastHisto2D(const FastHisto2D &other) : } FastHisto2D::T FastHisto2D::GetAt(const T &x, const T &y) const { - auto matchx = std::lower_bound(binEdgesX_.begin(), binEdgesX_.end(), x); + auto matchx = std::upper_bound(binEdgesX_.begin(), binEdgesX_.end(), x); if (matchx == binEdgesX_.begin() || matchx == binEdgesX_.end()) return T(0.0); int ix = (matchx - binEdgesX_.begin() - 1); - auto matchy = std::lower_bound(binEdgesY_.begin(), binEdgesY_.end(), y); + auto matchy = std::upper_bound(binEdgesY_.begin(), binEdgesY_.end(), y); if (matchy == binEdgesY_.begin() || matchy == binEdgesY_.end()) return T(0.0); int iy = (matchy - binEdgesY_.begin() - 1); return values_[ix * binY_ + iy]; @@ -186,7 +186,7 @@ FastHisto2D::T FastHisto2D::GetMaxOnXY() const { FastHisto2D::T FastHisto2D::GetMaxOnX(const T &y) const { - auto matchy = std::lower_bound(binEdgesY_.begin(), binEdgesY_.end(), y); + auto matchy = std::upper_bound(binEdgesY_.begin(), binEdgesY_.end(), y); if (matchy == binEdgesY_.begin() || matchy == binEdgesY_.end()) return T(0.0); int iy = (matchy - binEdgesY_.begin() - 1); T ret = 0.0; @@ -197,7 +197,7 @@ FastHisto2D::T FastHisto2D::GetMaxOnX(const T &y) const { } FastHisto2D::T FastHisto2D::GetMaxOnY(const T &x) const { - auto matchx = std::lower_bound(binEdgesX_.begin(), binEdgesX_.end(), x); + auto matchx = std::upper_bound(binEdgesX_.begin(), binEdgesX_.end(), x); if (matchx == binEdgesX_.begin() || matchx == binEdgesX_.end()) return T(0.0); int ix = (matchx - binEdgesX_.begin() - 1); return *std::max( &values_[ix * binY_], &values_[(ix+1) * binY_] ); @@ -247,13 +247,13 @@ FastHisto3D::FastHisto3D(const FastHisto3D &other) : } FastHisto3D::T FastHisto3D::GetAt(const T &x, const T &y, const T &z) const { - auto matchx = std::lower_bound(binEdgesX_.begin(), binEdgesX_.end(), x); + auto matchx = std::upper_bound(binEdgesX_.begin(), binEdgesX_.end(), x); if (matchx == binEdgesX_.begin() || matchx == binEdgesX_.end()) return T(0.0); int ix = (matchx - binEdgesX_.begin() - 1); - auto matchy = std::lower_bound(binEdgesY_.begin(), binEdgesY_.end(), y); + auto matchy = std::upper_bound(binEdgesY_.begin(), binEdgesY_.end(), y); if (matchy == binEdgesY_.begin() || matchy == binEdgesY_.end()) return T(0.0); int iy = (matchy - binEdgesY_.begin() - 1); - auto matchz = std::lower_bound(binEdgesZ_.begin(), binEdgesZ_.end(), z); + auto matchz = std::upper_bound(binEdgesZ_.begin(), binEdgesZ_.end(), z); if (matchz == binEdgesZ_.begin() || matchz == binEdgesZ_.end()) return T(0.0); int iz = (matchz - binEdgesZ_.begin() - 1); return values_[ix * binY_ *binZ_ +binZ_*iy + iz]; diff --git a/test/testCreateNLL.cxx b/test/testCreateNLL.cxx index 6c0488f85cc..6a1560aff65 100644 --- a/test/testCreateNLL.cxx +++ b/test/testCreateNLL.cxx @@ -31,213 +31,219 @@ namespace { -bool loadSnapshotIfExists(RooWorkspace &ws, const char *name) { - if (!name) - return false; - if (!ws.getSnapshot(name)) - return false; - return ws.loadSnapshot(name); -} + bool loadSnapshotIfExists(RooWorkspace& ws, const char* name) { + if (!name) + return false; + if (!ws.getSnapshot(name)) + return false; + return ws.loadSnapshot(name); + } -void configureCascadeMinimizerState(RooWorkspace &ws, RooStats::ModelConfig &mc) { - // Adapted from https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/master/src/Combine.cc#L1146-L1203 - auto &cfg = CascadeMinimizerGlobalConfigs::O(); + void configureCascadeMinimizerState(RooWorkspace& ws, RooStats::ModelConfig& mc) { + // Adapted from https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/master/src/Combine.cc#L1146-L1203 + auto& cfg = CascadeMinimizerGlobalConfigs::O(); - cfg.parametersOfInterest = RooArgList(); - if (const RooArgSet *poi = mc.GetParametersOfInterest()) { - for (RooAbsArg *arg : *poi) { - auto *var = dynamic_cast(arg); - if (!var) - continue; - cfg.parametersOfInterest.add(*var); + cfg.parametersOfInterest = RooArgList(); + if (const RooArgSet* poi = mc.GetParametersOfInterest()) { + for (RooAbsArg* arg : *poi) { + auto* var = dynamic_cast(arg); + if (!var) + continue; + cfg.parametersOfInterest.add(*var); + } + } + + cfg.nuisanceParameters = RooArgList(); + if (const RooArgSet* nuis = mc.GetNuisanceParameters()) { + for (RooAbsArg* arg : *nuis) { + auto* var = dynamic_cast(arg); + if (!var) + continue; + if (!var->isConstant()) + cfg.nuisanceParameters.add(*var); + } } - } - cfg.nuisanceParameters = RooArgList(); - if (const RooArgSet *nuis = mc.GetNuisanceParameters()) { - for (RooAbsArg *arg : *nuis) { - auto *var = dynamic_cast(arg); + cfg.allFloatingParameters = RooArgList(); + RooArgSet allVars(ws.allVars()); + for (RooAbsArg* arg : allVars) { + auto* var = dynamic_cast(arg); if (!var) continue; if (!var->isConstant()) - cfg.nuisanceParameters.add(*var); + cfg.allFloatingParameters.add(*var); } - } - cfg.allFloatingParameters = RooArgList(); - RooArgSet allVars(ws.allVars()); - for (RooAbsArg *arg : allVars) { - auto *var = dynamic_cast(arg); - if (!var) - continue; - if (!var->isConstant()) - cfg.allFloatingParameters.add(*var); - } + cfg.pdfCategories = RooArgList(); + cfg.allRooMultiPdfParams = RooArgList(); + cfg.allRooMultiPdfs = RooArgList(); - cfg.pdfCategories = RooArgList(); - cfg.allRooMultiPdfParams = RooArgList(); - cfg.allRooMultiPdfs = RooArgList(); - - if (auto *discreteParameters = dynamic_cast(ws.genobj("discreteParams")); - discreteParameters != nullptr) { - for (RooAbsArg *arg : *discreteParameters) { - auto *cat = dynamic_cast(arg); - if (!cat) - continue; - if (!cfg.pdfCategories.containsInstance(*cat)) - cfg.pdfCategories.add(*cat); - } - } else if (runtimedef::get("ADD_DISCRETE_FALLBACK")) { - RooArgSet categories(ws.allCats()); - for (RooAbsArg *arg : categories) { - auto *cat = dynamic_cast(arg); - if (!cat) - continue; - if (std::string(cat->GetName()).find("pdfindex") == std::string::npos) - continue; - if (!cfg.pdfCategories.containsInstance(*cat)) - cfg.pdfCategories.add(*cat); + if (auto* discreteParameters = dynamic_cast(ws.genobj("discreteParams")); + discreteParameters != nullptr) { + for (RooAbsArg* arg : *discreteParameters) { + auto* cat = dynamic_cast(arg); + if (!cat) + continue; + if (!cfg.pdfCategories.containsInstance(*cat)) + cfg.pdfCategories.add(*cat); + } + } else if (runtimedef::get("ADD_DISCRETE_FALLBACK")) { + RooArgSet categories(ws.allCats()); + for (RooAbsArg* arg : categories) { + auto* cat = dynamic_cast(arg); + if (!cat) + continue; + if (std::string(cat->GetName()).find("pdfindex") == std::string::npos) + continue; + if (!cfg.pdfCategories.containsInstance(*cat)) + cfg.pdfCategories.add(*cat); + } } - } - - if (cfg.pdfCategories.getSize() > 0) { - // Mirrors https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/master/src/Combine.cc#L1205-L1230 - RooArgSet clients; - utils::getClients(cfg.pdfCategories, ws.allPdfs(), clients); - for (RooAbsArg *arg : clients) { - auto *multi = dynamic_cast(arg); - if (!multi) - continue; - if (!cfg.allRooMultiPdfs.containsInstance(*multi)) - cfg.allRooMultiPdfs.add(*multi); - std::unique_ptr pdfParams{multi->getParameters((const RooArgSet *)nullptr)}; - for (RooAbsArg *a : *pdfParams) { - auto *var = dynamic_cast(a); - if (!var || var->isConstant()) + if (cfg.pdfCategories.getSize() > 0) { + // Mirrors https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/master/src/Combine.cc#L1205-L1230 + RooArgSet clients; + utils::getClients(cfg.pdfCategories, ws.allPdfs(), clients); + for (RooAbsArg* arg : clients) { + auto* multi = dynamic_cast(arg); + if (!multi) continue; - if (!cfg.allRooMultiPdfParams.containsInstance(*var)) - cfg.allRooMultiPdfParams.add(*var); + if (!cfg.allRooMultiPdfs.containsInstance(*multi)) + cfg.allRooMultiPdfs.add(*multi); + + std::unique_ptr pdfParams{multi->getParameters((const RooArgSet*)nullptr)}; + for (RooAbsArg* a : *pdfParams) { + auto* var = dynamic_cast(a); + if (!var || var->isConstant()) + continue; + if (!cfg.allRooMultiPdfParams.containsInstance(*var)) + cfg.allRooMultiPdfParams.add(*var); + } } } } -} -enum class BackendComparison { Cpu, Codegen }; - -struct FitSummary { - std::map poiValues; - std::map parameterErrors; - std::unique_ptr fitResult; - double minNll = std::numeric_limits::quiet_NaN(); -}; - -std::optional runFit(RooWorkspace &workspace, - RooStats::ModelConfig &modelConfig, - RooAbsPdf &pdf, - RooAbsData &data, - const RooArgSet &poiSet, - const std::string &backend, - std::string *error = nullptr) { - std:: cout << "Running fit with backend " << backend << std::endl; - if (!loadSnapshotIfExists(workspace, "clean")) { - for (RooAbsArg *arg : poiSet) { - if (auto *var = dynamic_cast(arg)) - var->setConstant(false); + enum class BackendComparison { Cpu, Codegen }; + + struct FitSummary { + std::map poiValues; + std::map parameterErrors; + std::unique_ptr fitResult; + double minNll = std::numeric_limits::quiet_NaN(); + }; + + std::optional runFit(RooWorkspace& workspace, + RooStats::ModelConfig& modelConfig, + RooAbsPdf& pdf, + RooAbsData& data, + const RooArgSet& poiSet, + const std::string& backend, + std::string* error = nullptr) { + std::cout << "Running fit with backend " << backend << std::endl; + if (!loadSnapshotIfExists(workspace, "clean")) { + for (RooAbsArg* arg : poiSet) { + if (auto* var = dynamic_cast(arg)) + var->setConstant(false); + } } - } - configureCascadeMinimizerState(workspace, modelConfig); - - const RooArgSet *nuisances = modelConfig.GetNuisanceParameters(); - std::unique_ptr constraintCopy; - const RooArgSet *constraintPtr = nullptr; - if (nuisances && nuisances->getSize() > 0) { - constraintCopy = std::make_unique(*nuisances); - constraintPtr = constraintCopy.get(); - } + configureCascadeMinimizerState(workspace, modelConfig); - Combine::setNllBackend(backend); - std::unique_ptr nll; - try { - nll = combineCreateNLL(pdf, data, constraintPtr, /*offset=*/true); - } catch (const std::exception &ex) { - if (error) *error = ex.what(); - return std::nullopt; - } - if (!nll) { - if (error) *error = "combineCreateNLL returned null pointer"; - return std::nullopt; - } + const RooArgSet* nuisances = modelConfig.GetNuisanceParameters(); + std::unique_ptr constraintCopy; + const RooArgSet* constraintPtr = nullptr; + if (nuisances && nuisances->getSize() > 0) { + constraintCopy = std::make_unique(*nuisances); + constraintPtr = constraintCopy.get(); + } - auto *primaryPoi = dynamic_cast(poiSet.first()); - CascadeMinimizer minim(*nll, CascadeMinimizer::Constrained, primaryPoi); - try { - if (!minim.minimize(/*verbose=*/0)) { - if (error) *error = "minimization failed"; + Combine::setNllBackend(backend); + std::unique_ptr nll; + try { + nll = combineCreateNLL(pdf, data, constraintPtr, /*offset=*/true); + if (backend == "codegen" || backend == "codegen_no_grad") { + RooFit::Experimental::writeCodegenDebugMacro(*nll, "codegen"); + } + } catch (const std::exception& ex) { + if (error) + *error = ex.what(); + return std::nullopt; + } + if (!nll) { + if (error) + *error = "combineCreateNLL returned null pointer"; return std::nullopt; } - } catch (const std::exception &ex) { - if (error) *error = ex.what(); - return std::nullopt; - } - minim.hesse(/*verbose=*/0); - auto fitResult = std::unique_ptr(minim.save()); + auto* primaryPoi = dynamic_cast(poiSet.first()); + CascadeMinimizer minim(*nll, CascadeMinimizer::Constrained, primaryPoi); + try { + if (!minim.minimize(/*verbose=*/0)) { + if (error) + *error = "minimization failed"; + return std::nullopt; + } + } catch (const std::exception& ex) { + if (error) + *error = ex.what(); + return std::nullopt; + } - FitSummary summary; - summary.minNll = nll->getVal(); - summary.fitResult = std::move(fitResult); + minim.hesse(/*verbose=*/0); + auto fitResult = std::unique_ptr(minim.save()); - for (RooAbsArg *arg : poiSet) { - auto *poiVar = dynamic_cast(arg); - if (!poiVar) - continue; - summary.poiValues[poiVar->GetName()] = poiVar->getVal(); - } + FitSummary summary; + summary.minNll = nll->getVal(); + summary.fitResult = std::move(fitResult); - if (summary.fitResult) { - const RooArgList &floats = summary.fitResult->floatParsFinal(); - for (int i = 0; i < floats.getSize(); ++i) { - auto *var = dynamic_cast(floats.at(i)); - if (!var) + for (RooAbsArg* arg : poiSet) { + auto* poiVar = dynamic_cast(arg); + if (!poiVar) continue; - const std::string name = var->GetName(); - const double err = var->getError(); - summary.parameterErrors[name] = err; - const double val = var->getVal(); - std::cout << " " << name << " = " << val; - if (!std::isnan(err)) { - std::cout << " +/-" << err; + summary.poiValues[poiVar->GetName()] = poiVar->getVal(); + } + + if (summary.fitResult) { + const RooArgList& floats = summary.fitResult->floatParsFinal(); + for (int i = 0; i < floats.getSize(); ++i) { + auto* var = dynamic_cast(floats.at(i)); + if (!var) + continue; + const std::string name = var->GetName(); + const double err = var->getError(); + summary.parameterErrors[name] = err; + const double val = var->getVal(); + std::cout << " " << name << " = " << val; + if (!std::isnan(err)) { + std::cout << " +/-" << err; + } + std::cout << '\n'; } - std::cout << '\n'; } + + return summary; } - return summary; -} + void compareFits(const FitSummary& baseline, + const FitSummary& candidate, + const std::string& label, + double relTolerance) { + for (const auto& entry : baseline.parameterErrors) { + const std::string& name = entry.first; + if (name.rfind("prop_", 0) != 0) + continue; + const auto it = candidate.parameterErrors.find(name); + ASSERT_TRUE(it != candidate.parameterErrors.end()) + << "Parameter '" << name << "' missing in backend '" << label << "'"; + + const double err1 = entry.second; + const double err2 = it->second; + const double diff = std::abs(err1 - err2); + const double scale = std::max({1.0, std::abs(err1), std::abs(err2)}); -void compareFits(const FitSummary &baseline, - const FitSummary &candidate, - const std::string &label, - double relTolerance) { - for (const auto &entry : baseline.parameterErrors) { - const std::string &name = entry.first; - if (name.rfind("prop_", 0) != 0) - continue; - const auto it = candidate.parameterErrors.find(name); - ASSERT_TRUE(it != candidate.parameterErrors.end()) - << "Parameter '" << name << "' missing in backend '" << label << "'"; - - const double err1 = entry.second; - const double err2 = it->second; - const double diff = std::abs(err1 - err2); - const double scale = std::max({1.0, std::abs(err1), std::abs(err2)}); - - EXPECT_TRUE(diff / scale <= relTolerance) - << "Parameter '" << name << "' error mismatch: combine=" << err1 - << ", " << label << "=" << err2; + EXPECT_TRUE(diff / scale <= relTolerance) + << "Parameter '" << name << "' error mismatch: combine=" << err1 << ", " << label << "=" << err2; + } } -} } // namespace @@ -255,142 +261,138 @@ class CreateNLLTest : public ::testing::TestWithParam { runtimedef::set("MINIMIZER_no_analytic", enableAnalytic ? 0 : 1); } - void TearDown() override { - runtimedef::set("MINIMIZER_no_analytic", oldAnalyticFlag_); - } + void TearDown() override { runtimedef::set("MINIMIZER_no_analytic", oldAnalyticFlag_); } void runComparison() { - CascadeMinimizer::initOptions(); - // Matches https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/master/bin/combine.cpp#L300-L317 - runtimedef::set("OPTIMIZE_BOUNDS", 1); - runtimedef::set("ADDNLL_RECURSIVE", 1); - runtimedef::set("ADDNLL_GAUSSNLL", 1); - runtimedef::set("ADDNLL_HISTNLL", 1); - runtimedef::set("ADDNLL_CBNLL", 1); - runtimedef::set("TMCSO_AdaptivePseudoAsimov", 1); - runtimedef::set("MINIMIZER_optimizeConst", 2); - runtimedef::set("MINIMIZER_rooFitOffset", 1); - runtimedef::set("ADDNLL_ROOREALSUM_FACTOR", 1); - runtimedef::set("ADDNLL_ROOREALSUM_NONORM", 1); - runtimedef::set("ADDNLL_ROOREALSUM_BASICINT", 1); - runtimedef::set("ADDNLL_ROOREALSUM_KEEPZEROS", 1); - runtimedef::set("ADDNLL_PRODNLL", 1); - runtimedef::set("ADDNLL_HFNLL", 1); - runtimedef::set("ADDNLL_HISTFUNCNLL", 1); - runtimedef::set("ADDNLL_ROOREALSUM_CHEAPPROD", 1); - - const std::string fileName = "template-analysis_shape_autoMCStats.root"; - const std::string workspaceName = "w"; - const std::string modelConfigName = "ModelConfig"; - const std::string dataName = "data_obs"; - - std::unique_ptr file(TFile::Open(fileName.c_str(), "READ")); - EXPECT_FALSE(!file || file->IsZombie()) << "ERROR: failed to open input file '" << fileName << "'.\n"; - - auto *workspace = dynamic_cast(file->Get(workspaceName.c_str())); - if (!workspace) { - std::cerr << "ERROR: workspace '" << workspaceName << "' not found in '" << fileName << "'.\n"; - file->ls(); - //return 2; - } + CascadeMinimizer::initOptions(); + // Matches https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/master/bin/combine.cpp#L300-L317 + runtimedef::set("OPTIMIZE_BOUNDS", 1); + runtimedef::set("ADDNLL_RECURSIVE", 1); + runtimedef::set("ADDNLL_GAUSSNLL", 1); + runtimedef::set("ADDNLL_HISTNLL", 1); + runtimedef::set("ADDNLL_CBNLL", 1); + runtimedef::set("TMCSO_AdaptivePseudoAsimov", 1); + runtimedef::set("MINIMIZER_optimizeConst", 2); + runtimedef::set("MINIMIZER_rooFitOffset", 1); + runtimedef::set("ADDNLL_ROOREALSUM_FACTOR", 1); + runtimedef::set("ADDNLL_ROOREALSUM_NONORM", 1); + runtimedef::set("ADDNLL_ROOREALSUM_BASICINT", 1); + runtimedef::set("ADDNLL_ROOREALSUM_KEEPZEROS", 1); + runtimedef::set("ADDNLL_PRODNLL", 1); + runtimedef::set("ADDNLL_HFNLL", 1); + runtimedef::set("ADDNLL_HISTFUNCNLL", 1); + runtimedef::set("ADDNLL_ROOREALSUM_CHEAPPROD", 1); + + const std::string fileName = "template-analysis_shape_autoMCStats.root"; + const std::string workspaceName = "w"; + const std::string modelConfigName = "ModelConfig"; + const std::string dataName = "data_obs"; + + std::unique_ptr file(TFile::Open(fileName.c_str(), "READ")); + EXPECT_FALSE(!file || file->IsZombie()) << "ERROR: failed to open input file '" << fileName << "'.\n"; + + auto* workspace = dynamic_cast(file->Get(workspaceName.c_str())); + if (!workspace) { + std::cerr << "ERROR: workspace '" << workspaceName << "' not found in '" << fileName << "'.\n"; + file->ls(); + //return 2; + } - // Matches https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/master/src/Combine.cc#L505-L603 - utils::check_inf_parameters(workspace->allVars(), /*verbosity=*/0); - loadSnapshotIfExists(*workspace, "clean"); + // Matches https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/master/src/Combine.cc#L505-L603 + utils::check_inf_parameters(workspace->allVars(), /*verbosity=*/0); + loadSnapshotIfExists(*workspace, "clean"); - auto *modelConfig = - dynamic_cast(workspace->genobj(modelConfigName.c_str())); - if (!modelConfig) { - std::cerr << "ERROR: ModelConfig '" << modelConfigName << "' not found in workspace '" - << workspaceName << "'.\n"; - //return 2; - } + auto* modelConfig = dynamic_cast(workspace->genobj(modelConfigName.c_str())); + if (!modelConfig) { + std::cerr << "ERROR: ModelConfig '" << modelConfigName << "' not found in workspace '" << workspaceName << "'.\n"; + //return 2; + } - RooAbsPdf *pdf = modelConfig->GetPdf(); - if (!pdf) { - std::cerr << "ERROR: ModelConfig '" << modelConfigName << "' does not define a pdf.\n"; - //return 2; - } + RooAbsPdf* pdf = modelConfig->GetPdf(); + if (!pdf) { + std::cerr << "ERROR: ModelConfig '" << modelConfigName << "' does not define a pdf.\n"; + //return 2; + } - RooAbsData *data = workspace->data(dataName.c_str()); - if (!data) { - std::cerr << "ERROR: dataset '" << dataName << "' not found in workspace '" << workspaceName - << "'.\n"; - //return 2; - } + RooAbsData* data = workspace->data(dataName.c_str()); + if (!data) { + std::cerr << "ERROR: dataset '" << dataName << "' not found in workspace '" << workspaceName << "'.\n"; + //return 2; + } - RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); - RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); + RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); + RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CountErrors); - const RooArgSet *poi = modelConfig->GetParametersOfInterest(); - if (!poi || poi->getSize() == 0) { - std::cerr << "ERROR: ModelConfig '" << modelConfigName << "' has no parameters of interest.\n"; - //return 2; - } - for (RooAbsArg *arg : *poi) { - auto *var = dynamic_cast(arg); - if (var) - var->setConstant(false); - } + const RooArgSet* poi = modelConfig->GetParametersOfInterest(); + if (!poi || poi->getSize() == 0) { + std::cerr << "ERROR: ModelConfig '" << modelConfigName << "' has no parameters of interest.\n"; + //return 2; + } + for (RooAbsArg* arg : *poi) { + auto* var = dynamic_cast(arg); + if (var) + var->setConstant(false); + } - const TestConfig cfg = GetParam(); + const TestConfig cfg = GetParam(); - std::string errMessage; - auto combineFit = runFit(*workspace, *modelConfig, *pdf, *data, *poi, "combine", &errMessage); - ASSERT_TRUE(combineFit) << "Failed to run fit with 'combine' backend: " << errMessage; + std::string errMessage; + auto combineFit = runFit(*workspace, *modelConfig, *pdf, *data, *poi, "combine", &errMessage); + ASSERT_TRUE(combineFit) << "Failed to run fit with 'combine' backend: " << errMessage; - errMessage.clear(); - const double relTolerance = 0.05; // 5% + errMessage.clear(); + const double relTolerance = 0.05; // 5% - if (cfg.backend == BackendComparison::Cpu) { - auto cpuFit = runFit(*workspace, *modelConfig, *pdf, *data, *poi, "cpu", &errMessage); - ASSERT_TRUE(cpuFit) << "Failed to run fit with 'cpu' backend: " << errMessage; - if (cfg.expectFailure) { - ADD_FAILURE() << "Expected compare failure for cpu backend but fit succeeded"; - } - compareFits(*combineFit, *cpuFit, "cpu", relTolerance); - } else { - auto codegenFit = runFit(*workspace, *modelConfig, *pdf, *data, *poi, "codegen", &errMessage); - if (!codegenFit) { + if (cfg.backend == BackendComparison::Cpu) { + auto cpuFit = runFit(*workspace, *modelConfig, *pdf, *data, *poi, "cpu", &errMessage); + ASSERT_TRUE(cpuFit) << "Failed to run fit with 'cpu' backend: " << errMessage; if (cfg.expectFailure) { - GTEST_SKIP() << "Expected failure for codegen backend: " << errMessage; - } else { - ADD_FAILURE() << "Failed to run fit with 'codegen' backend: " << errMessage; + ADD_FAILURE() << "Expected compare failure for cpu backend but fit succeeded"; } - return; - } - if (cfg.expectFailure) { - ADD_FAILURE() << "Expected compare failure for codegen backend but fit succeeded"; + compareFits(*combineFit, *cpuFit, "cpu", relTolerance); + } else { + // we don't use the analytic gradient yet. For now we just want to + // validate the generated code. + auto codegenFit = runFit(*workspace, *modelConfig, *pdf, *data, *poi, "codegen_no_grad", &errMessage); + if (!codegenFit) { + if (cfg.expectFailure) { + GTEST_SKIP() << "Expected failure for codegen backend: " << errMessage; + } else { + ADD_FAILURE() << "Failed to run fit with 'codegen' backend: " << errMessage; + } + return; + } + if (cfg.expectFailure) { + ADD_FAILURE() << "Expected compare failure for codegen backend but fit succeeded"; + } + compareFits(*combineFit, *codegenFit, "codegen", relTolerance); } - compareFits(*combineFit, *codegenFit, "codegen", relTolerance); - } } private: int oldAnalyticFlag_ = 0; }; -TEST_P(CreateNLLTest, BackendConsistency) { - runComparison(); -} +TEST_P(CreateNLLTest, BackendConsistency) { runComparison(); } -static std::string AnalyticParamName(const testing::TestParamInfo &info) { - const char *backend = info.param.backend == BackendComparison::Cpu ? "Cpu" : "Codegen"; - const char *analytic = info.param.analyticBarlowBeeston ? "AnalyticOn" : "AnalyticOff"; - const char *expect = info.param.expectFailure ? "ExpectFail" : "ExpectPass"; +static std::string AnalyticParamName(const testing::TestParamInfo& info) { + const char* backend = info.param.backend == BackendComparison::Cpu ? "Cpu" : "Codegen"; + const char* analytic = info.param.analyticBarlowBeeston ? "AnalyticOn" : "AnalyticOff"; + const char* expect = info.param.expectFailure ? "ExpectFail" : "ExpectPass"; return std::string(backend) + "_" + analytic + "_" + expect; } INSTANTIATE_TEST_SUITE_P( - // Alternative evaluation backends are only supported from ROOT 6.30. +// Alternative evaluation backends are only supported from ROOT 6.30. #if ROOT_VERSION_CODE < ROOT_VERSION(6, 30, 0) DISABLED_BarlowBeestonMatrix, #else BarlowBeestonMatrix, #endif CreateNLLTest, - ::testing::Values(TestConfig{BackendComparison::Cpu, true, false}, - TestConfig{BackendComparison::Codegen, true, true}, - TestConfig{BackendComparison::Cpu, false, false}, - TestConfig{BackendComparison::Codegen, false, true}), + ::testing::Values( + //TestConfig{BackendComparison::Cpu, true, false}, + //TestConfig{BackendComparison::Codegen, true, true}, + TestConfig{BackendComparison::Cpu, false, false}, + TestConfig{BackendComparison::Codegen, false, true}), AnalyticParamName);