From b33194cad30619d6dde3618d313e1fd4e4e46142 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sat, 8 Jun 2024 13:20:48 +0200 Subject: [PATCH] First RooFit AD integration These are the remaining developments from the RooFit ICHEP 2024 AD branch. --- interface/AsymPow.h | 5 + interface/CombineCodegenImpl.h | 18 +++ interface/HZZ4LRooPdfs.h | 10 ++ interface/ProcessNormalization.h | 6 + interface/VerticalInterpHistPdf.h | 9 +- scripts/fitRooFitAD.py | 38 ++++++ src/CombineCodegenImpl.cxx | 213 ++++++++++++++++++++++++++++++ src/HZZ4LRooPdfs.cc | 21 ++- src/VerticalInterpHistPdf.cc | 2 + 9 files changed, 320 insertions(+), 2 deletions(-) create mode 100644 interface/CombineCodegenImpl.h create mode 100644 scripts/fitRooFitAD.py create mode 100644 src/CombineCodegenImpl.cxx diff --git a/interface/AsymPow.h b/interface/AsymPow.h index 6cda05c5522..552c9829941 100644 --- a/interface/AsymPow.h +++ b/interface/AsymPow.h @@ -4,6 +4,7 @@ #include #include +#include "CombineCodegenImpl.h" //_________________________________________________ /* @@ -28,6 +29,8 @@ class AsymPow : public RooAbsReal { TObject * clone(const char *newname) const override { return new AsymPow(*this, newname); } + COMBINE_DECLARE_TRANSLATE; + RooAbsReal const &kappaLow() const { return kappaLow_.arg(); } RooAbsReal const &kappaHigh() const { return kappaHigh_.arg(); } RooAbsReal const &theta() const { return theta_.arg(); } @@ -42,4 +45,6 @@ class AsymPow : public RooAbsReal { ClassDefOverride(AsymPow,1) // Asymmetric power }; +COMBINE_DECLARE_CODEGEN_IMPL(AsymPow); + #endif diff --git a/interface/CombineCodegenImpl.h b/interface/CombineCodegenImpl.h new file mode 100644 index 00000000000..d1795eb1520 --- /dev/null +++ b/interface/CombineCodegenImpl.h @@ -0,0 +1,18 @@ +#ifndef HiggsAnalysis_CombinedLimit_CombineCodegenImpl_h +#define HiggsAnalysis_CombinedLimit_CombineCodegenImpl_h + +#include // for ROOT_VERSION + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,35,0) +# define COMBINE_DECLARE_CODEGEN_IMPL(CLASS_NAME) \ + namespace RooFit { namespace Experimental { void codegenImpl(CLASS_NAME &arg, CodegenContext &ctx); }} +# define COMBINE_DECLARE_TRANSLATE +#elif ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) +# define COMBINE_DECLARE_CODEGEN_IMPL(CLASS_NAME) +# define COMBINE_DECLARE_TRANSLATE void translate(RooFit::Detail::CodeSquashContext &ctx) const override; +#else +# define COMBINE_DECLARE_CODEGEN_IMPL(_) +# define COMBINE_DECLARE_TRANSLATE +#endif + +#endif diff --git a/interface/HZZ4LRooPdfs.h b/interface/HZZ4LRooPdfs.h index df1f7e2e3ae..fbae926a4d8 100644 --- a/interface/HZZ4LRooPdfs.h +++ b/interface/HZZ4LRooPdfs.h @@ -10,6 +10,8 @@ #include "TDirectory.h" #include "TH2F.h" +#include // for ROOT_VERSION + #include #include #include @@ -136,6 +138,10 @@ class RooqqZZPdf_v2 : public RooAbsPdf { RooqqZZPdf_v2(const RooqqZZPdf_v2& other, const char* name=0) ; TObject* clone(const char* newname) const override { return new RooqqZZPdf_v2(*this,newname); } inline ~RooqqZZPdf_v2() override { } + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) + std::unique_ptr compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const override; +#endif protected: @@ -299,6 +305,10 @@ class RooggZZPdf_v2 : public RooAbsPdf { RooggZZPdf_v2(const RooggZZPdf_v2& other, const char* name=0) ; TObject* clone(const char* newname) const override { return new RooggZZPdf_v2(*this,newname); } inline ~RooggZZPdf_v2() override { } + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) + std::unique_ptr compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const override; +#endif protected: diff --git a/interface/ProcessNormalization.h b/interface/ProcessNormalization.h index edf2ec3f15c..a3444c5d2e7 100644 --- a/interface/ProcessNormalization.h +++ b/interface/ProcessNormalization.h @@ -4,6 +4,8 @@ #include #include "RooListProxy.h" +#include "CombineCodegenImpl.h" + //_________________________________________________ /* BEGIN_HTML @@ -28,6 +30,8 @@ class ProcessNormalization : public RooAbsReal { void addOtherFactor(RooAbsReal &factor) ; void dump() const ; + COMBINE_DECLARE_TRANSLATE; + double nominalValue() const { return nominalValue_; } std::vector const &logKappa() const { return logKappa_; } RooArgList const &thetaList() const { return thetaList_; } @@ -57,4 +61,6 @@ class ProcessNormalization : public RooAbsReal { ClassDefOverride(ProcessNormalization,1) // Process normalization interpolator }; +COMBINE_DECLARE_CODEGEN_IMPL(ProcessNormalization); + #endif diff --git a/interface/VerticalInterpHistPdf.h b/interface/VerticalInterpHistPdf.h index c90fdc051ff..ccfc34aad26 100644 --- a/interface/VerticalInterpHistPdf.h +++ b/interface/VerticalInterpHistPdf.h @@ -10,7 +10,8 @@ #include "TH1.h" #include "SimpleCacheSentry.h" #include "FastTemplate_Old.h" -#include + +#include "CombineCodegenImpl.h" class FastVerticalInterpHistPdf; class FastVerticalInterpHistPdf2Base; @@ -328,6 +329,8 @@ class FastVerticalInterpHistPdf2 : public FastVerticalInterpHistPdf2Base { FastHisto const& cache() const { return _cache; } + COMBINE_DECLARE_TRANSLATE; + FastHisto const &cacheNominal() const { return _cacheNominal; } friend class FastVerticalInterpHistPdf2V; @@ -387,6 +390,8 @@ class FastVerticalInterpHistPdf2D2 : public FastVerticalInterpHistPdf2Base { Double_t evaluate() const override ; + COMBINE_DECLARE_TRANSLATE; + FastHisto2D const &cacheNominal() const { return _cacheNominal; } protected: @@ -459,5 +464,7 @@ class FastVerticalInterpHistPdf3D : public FastVerticalInterpHistPdfBase { }; +COMBINE_DECLARE_CODEGEN_IMPL(FastVerticalInterpHistPdf2); +COMBINE_DECLARE_CODEGEN_IMPL(FastVerticalInterpHistPdf2D2); #endif diff --git a/scripts/fitRooFitAD.py b/scripts/fitRooFitAD.py new file mode 100644 index 00000000000..0acf47366c9 --- /dev/null +++ b/scripts/fitRooFitAD.py @@ -0,0 +1,38 @@ +import numpy as np +import ROOT +import argparse + +parser = argparse.ArgumentParser() +parser.add_argument("--input", "-i", help="input ws file") +parser.add_argument("--backend", help='set evaluation backend, default is "combine"') + +args = parser.parse_args() + +with ROOT.TFile.Open(args.input) as f: + ws = f.Get("w") + +global_observables = ws.set("globalObservables") +constrain = ws.set("nuisances") + +pdf = ws["model_s"] +data = ws["data_obs"] + +# To use the evaluation backends of RooFit +ROOT.gInterpreter.Declare("#include ") + +ROOT.Combine.nllBackend_ = args.backend + +# Change verbosity +ROOT.RooMsgService.instance().getStream(1).removeTopic(ROOT.RooFit.Minimization) +ROOT.RooMsgService.instance().getStream(1).removeTopic(ROOT.RooFit.Fitting) + +ROOT.gInterpreter.Declare("#include ") + +nll = pdf.createNLL(data, Constrain=constrain, GlobalObservables=global_observables, Offset="initial") + +cfg = ROOT.RooMinimizer.Config() +minim = ROOT.RooMinimizer(nll, cfg) +minim.setEps(1.0) +minim.setStrategy(0) +minim.minimize("Minuit2", "") +minim.save().Print() diff --git a/src/CombineCodegenImpl.cxx b/src/CombineCodegenImpl.cxx new file mode 100644 index 00000000000..405a9695003 --- /dev/null +++ b/src/CombineCodegenImpl.cxx @@ -0,0 +1,213 @@ +#include "../interface/CombineCodegenImpl.h" + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) + +#include "../interface/AsymPow.h" +#include "../interface/ProcessNormalization.h" +#include "../interface/VerticalInterpHistPdf.h" + +#include + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,35,0) +namespace RooFit { +namespace Experimental { +# define CODEGEN_IMPL(CLASS_NAME) void codegenImpl(CLASS_NAME &arg0, CodegenContext &ctx) +# define ARG_VAR auto &arg = arg0; +#else +# define CODEGEN_IMPL(CLASS_NAME) void CLASS_NAME::translate(RooFit::Detail::CodeSquashContext &ctx) const +# define ARG_VAR auto &arg = *this; +#endif + + +CODEGEN_IMPL(AsymPow) { + ARG_VAR; + ctx.addResult(&arg, + ctx.buildCall("RooFit::Detail::MathFuncs::asymPow", arg.theta(), arg.kappaLow(), arg.kappaHigh())); +} + +CODEGEN_IMPL(ProcessNormalization) { + ARG_VAR; + + std::vector logAsymmKappaLow; + std::vector logAsymmKappaHigh; + logAsymmKappaLow.reserve(arg.logAsymmKappa().size()); + logAsymmKappaHigh.reserve(arg.logAsymmKappa().size()); + for (auto [lo, hi] : arg.logAsymmKappa()) { + logAsymmKappaLow.push_back(lo); + logAsymmKappaHigh.push_back(hi); + } + + ctx.addResult(&arg, + ctx.buildCall("RooFit::Detail::MathFuncs::processNormalization", + arg.nominalValue(), + arg.thetaList().size(), + arg.asymmThetaList().size(), + arg.otherFactorList().size(), + arg.thetaList(), + arg.logKappa(), + arg.asymmThetaList(), + logAsymmKappaLow, + logAsymmKappaHigh, + arg.otherFactorList())); +} + +CODEGEN_IMPL(FastVerticalInterpHistPdf2) { + ARG_VAR; + + if (arg.smoothAlgo() < 0) { + throw std::runtime_error("We only support smoothAlgo >= 0"); + } + + RooRealVar const &xVar = arg.x(); + + int numBins = xVar.numBins(); + + std::vector nominalVec(numBins); + std::vector widthVec(numBins); + std::vector morphsVecSum; + std::vector morphsVecDiff; + + auto const &cacheNominal = arg.cacheNominal(); + + for (int i = 0; i < numBins; ++i) { + nominalVec[i] = cacheNominal.GetBinContent(i); + widthVec[i] = cacheNominal.GetWidth(i); + } + + std::size_t nCoefs = arg.coefList().size(); + + morphsVecSum.reserve(numBins * nCoefs); + morphsVecDiff.reserve(numBins * nCoefs); + auto const &morphs = arg.morphs(); + for (unsigned int j = 0; j < nCoefs; ++j) { + for (int i = 0; i < numBins; ++i) { + morphsVecSum.push_back(morphs[j].sum[i]); + morphsVecDiff.push_back(morphs[j].diff[i]); + } + } + + for (int i = 0; i < numBins; ++i) { + nominalVec[i] = cacheNominal.GetBinContent(i); + } + + // The bin index part + // We also have to assert that x is uniformely binned! + if (!dynamic_cast(&xVar.getBinning())) { + throw std::runtime_error("We only support uniform binning!"); + } + double xLow = xVar.getMin(); + double xHigh = xVar.getMax(); + std::string binIdx = ctx.buildCall("RooFit::Detail::MathFuncs::getUniformBinning", xLow, xHigh, xVar, numBins); + + std::string arrName = ctx.getTmpVarName(); + + std::stringstream code; + code << "double " << arrName << "[" << numBins << "];\n"; + code << ctx.buildCall("RooFit::Detail::MathFuncs::fastVerticalInterpHistPdf2", + numBins, + nCoefs, + arg.coefList(), + nominalVec, + widthVec, + morphsVecSum, + morphsVecDiff, + arg.smoothRegion(), + arrName) + + ";\n"; + + ctx.addToCodeBody(code.str(), true); + ctx.addResult(&arg, arrName + "[" + binIdx + "]"); +} + +CODEGEN_IMPL(FastVerticalInterpHistPdf2D2) { + ARG_VAR; + + if (arg.smoothAlgo() < 0) { + throw std::runtime_error("We only support smoothAlgo >= 0"); + } + + if (!arg.conditional()) { + throw std::runtime_error("We only support conditional == true"); + } + + RooRealVar const &xVar = arg.x(); + RooRealVar const &yVar = arg.y(); + + // We also have to assert that x and y are uniformely binned! + if (!dynamic_cast(&xVar.getBinning())) { + throw std::runtime_error("We only support uniform binning!"); + } + if (!dynamic_cast(&yVar.getBinning())) { + throw std::runtime_error("We only support uniform binning!"); + } + + auto const &cacheNominal = arg.cacheNominal(); + + int numBinsX = cacheNominal.binX(); + int numBinsY = cacheNominal.binY(); + int numBins = numBinsY * numBinsY; + + std::vector nominalVec(numBins); + std::vector widthVec(numBins); + std::vector morphsVecSum; + std::vector morphsVecDiff; + + for (int i = 0; i < numBins; ++i) { + nominalVec[i] = cacheNominal.GetBinContent(i); + widthVec[i] = cacheNominal.GetWidth(i); + } + + std::size_t nCoefs = arg.coefList().size(); + + morphsVecSum.reserve(numBins * nCoefs); + morphsVecDiff.reserve(numBins * nCoefs); + auto const &morphs = arg.morphs(); + for (unsigned int j = 0; j < nCoefs; ++j) { + for (int i = 0; i < numBins; ++i) { + morphsVecSum.push_back(morphs[j].sum[i]); + morphsVecDiff.push_back(morphs[j].diff[i]); + } + } + + for (int i = 0; i < numBins; ++i) { + nominalVec[i] = cacheNominal.GetBinContent(i); + } + + // The bin index part + double xLow = xVar.getMin(); + double xHigh = xVar.getMax(); + std::string binIdxX = ctx.buildCall("RooFit::Detail::MathFuncs::getUniformBinning", xLow, xHigh, arg.x(), numBinsX); + double yLow = yVar.getMin(); + double yHigh = yVar.getMax(); + std::string binIdxY = ctx.buildCall("RooFit::Detail::MathFuncs::getUniformBinning", yLow, yHigh, arg.y(), numBinsY); + + std::stringstream binIdx; + binIdx << "(" << binIdxY << " + " << yVar.numBins() << " * " << binIdxX << ")"; + + std::string arrName = ctx.getTmpVarName(); + + std::stringstream code; + code << "double " << arrName << "[" << (numBinsX * numBinsY) << "];\n"; + code << ctx.buildCall("RooFit::Detail::MathFuncs::fastVerticalInterpHistPdf2D2", + numBinsX, + numBinsY, + nCoefs, + arg.coefList(), + nominalVec, + widthVec, + morphsVecSum, + morphsVecDiff, + arg.smoothRegion(), + arrName) + + ";\n"; + + ctx.addToCodeBody(code.str(), true); + ctx.addResult(&arg, arrName + "[" + binIdx.str() + "]"); +} + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,35,0) +} // namespace RooFit +} // namespace Experimental +#endif + +#endif // ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) diff --git a/src/HZZ4LRooPdfs.cc b/src/HZZ4LRooPdfs.cc index 613031da83b..302638513b3 100644 --- a/src/HZZ4LRooPdfs.cc +++ b/src/HZZ4LRooPdfs.cc @@ -4,7 +4,7 @@ #include "RooAbsReal.h" #include "RooRealVar.h" #include "RooVoigtian.h" -//#include "RooComplex.h" +#include "RooGenericPdf.h" #include "RooMath.h" #include "RooAbsCategory.h" #include @@ -2126,6 +2126,15 @@ Double_t RooqqZZPdf_v2::evaluate() const } +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) +std::unique_ptr RooqqZZPdf_v2::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const +{ + return RooGenericPdf{GetName(), + "(.5+.5*TMath::Erf((x[0]-x[1])/x[2]))*(x[4]/(1+exp((x[0]-x[1])/x[3])))+(.5+.5*TMath::Erf((x[0]-x[5])/x[6]))*(x[8]/(1+exp((x[0]-x[5])/x[7]))+x[10]/(1+exp((x[0]-x[5])/x[9])))+(.5+.5*TMath::Erf((x[0]-x[11])/x[12]))*(x[14]/(1+exp((x[0]-x[11])/x[13])) )", + {*m4l, *a0, *a1, *a2, *a3, *a4, *a5, *a6, *a7, *a8, *a9, *a10, *a11, *a12, *a13} + }.compileForNormSet(normSet, ctx); +} +#endif /************RooggZZPdf_v2***********/ @@ -2189,6 +2198,16 @@ Double_t RooggZZPdf_v2::evaluate() const return ZZ; } +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) +std::unique_ptr RooggZZPdf_v2::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const +{ + return RooGenericPdf{GetName(), + "(.5+.5*TMath::Erf((x[0]-x[1])/x[2]))*(x[4]/(1+exp((x[0]-x[1])/x[3]))) + (.5+.5*TMath::Erf((x[0]-x[5])/x[6]))*(x[8]/(1+exp((x[0]-x[5])/x[7]))+x[10]/(1+exp((x[0]-x[5])/x[9])))", + {*m4l, *a0, *a1, *a2, *a3, *a4, *a5, *a6, *a7, *a8, *a9} + }.compileForNormSet(normSet, ctx); +} +#endif + //// ------- RooVBFZZPdf ------ diff --git a/src/VerticalInterpHistPdf.cc b/src/VerticalInterpHistPdf.cc index ff9ef2ec65a..df2ac8ee5ed 100644 --- a/src/VerticalInterpHistPdf.cc +++ b/src/VerticalInterpHistPdf.cc @@ -10,6 +10,8 @@ #include "RooMsgService.h" #include "RooAbsData.h" +#include "RooUniformBinning.h" + //#define TRACE_CALLS #ifdef TRACE_CALLS #include "../interface/ProfilingTools.h"