diff --git a/roofit/roofit/inc/RooExponential.h b/roofit/roofit/inc/RooExponential.h index c9e646d02315b1c11930de483d6c2a9c388e03f4..a8e693cfc81cc1249752c56d0aba56cf0774f171 100644 --- a/roofit/roofit/inc/RooExponential.h +++ b/roofit/roofit/inc/RooExponential.h @@ -34,6 +34,10 @@ public: Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=nullptr) const override; double analyticalIntegral(Int_t code, const char* rangeName=nullptr) const override; + void translate(RooFit::Detail::CodeSquashContext &ctx) const override; + std::string + buildCallToAnalyticIntegral(Int_t code, const char *rangeName, RooFit::Detail::CodeSquashContext &ctx) const override; + protected: RooRealProxy x; RooRealProxy c; diff --git a/roofit/roofit/src/RooExponential.cxx b/roofit/roofit/src/RooExponential.cxx index a4eef6c6616b76b3a724981b61a8f0a8ba38f4ee..1311b3cf329725142d334adf2c269d6bf71bcdc8 100644 --- a/roofit/roofit/src/RooExponential.cxx +++ b/roofit/roofit/src/RooExponential.cxx @@ -30,11 +30,10 @@ range and values of the arguments. #include "RooRealVar.h" #include "RooBatchCompute.h" +#include <RooFit/Detail/AnalyticalIntegrals.h> #include <cmath> -using namespace std; - ClassImp(RooExponential); //////////////////////////////////////////////////////////////////////////////// @@ -57,7 +56,7 @@ RooExponential::RooExponential(const RooExponential& other, const char* name) : //////////////////////////////////////////////////////////////////////////////// double RooExponential::evaluate() const{ - return exp(c*x); + return std::exp(c*x); } //////////////////////////////////////////////////////////////////////////////// @@ -85,10 +84,25 @@ double RooExponential::analyticalIntegral(Int_t code, const char* rangeName) con auto& constant = code == 1 ? c : x; auto& integrand = code == 1 ? x : c; - if (constant == 0.0) { - return integrand.max(rangeName) - integrand.min(rangeName); - } + return RooFit::Detail::AnalyticalIntegrals::exponentialIntegral(integrand.min(rangeName), integrand.max(rangeName), constant); +} + +//////////////////////////////////////////////////////////////////////////////// + +void RooExponential::translate(RooFit::Detail::CodeSquashContext &ctx) const +{ + // Build a call to the stateless gaussian defined later. + ctx.addResult(this, "std::exp(" + ctx.getResult(c) + " * " + ctx.getResult(x) + ")"); +} + +//////////////////////////////////////////////////////////////////////////////// + +std::string RooExponential::buildCallToAnalyticIntegral(Int_t code, const char *rangeName, + RooFit::Detail::CodeSquashContext &ctx) const +{ + auto& constant = code == 1 ? c : x; + auto& integrand = code == 1 ? x : c; - return (exp(constant*integrand.max(rangeName)) - exp(constant*integrand.min(rangeName))) - / constant; + return ctx.buildCall("RooFit::Detail::AnalyticalIntegrals::exponentialIntegral", + integrand.min(rangeName), integrand.max(rangeName), constant); } diff --git a/roofit/roofit/src/RooGaussian.cxx b/roofit/roofit/src/RooGaussian.cxx index 3cc5e12b358c2d2d2793a59c8053891e89ed6b38..7083b7731a9395f31ed3c1e325a57f60acf8c0b6 100644 --- a/roofit/roofit/src/RooGaussian.cxx +++ b/roofit/roofit/src/RooGaussian.cxx @@ -24,7 +24,6 @@ Plain Gaussian p.d.f #include "RooBatchCompute.h" #include "RooHelpers.h" #include "RooRandom.h" -#include "RooNumber.h" #include <RooFit/Detail/AnalyticalIntegrals.h> #include <RooFit/Detail/EvaluateFuncs.h> @@ -58,7 +57,7 @@ RooGaussian::RooGaussian(const RooGaussian& other, const char* name) : double RooGaussian::evaluate() const { - return RooFit::Detail::EvaluateFuncs::gaussEvaluate(x, mean, sigma); + return RooFit::Detail::EvaluateFuncs::gaussianEvaluate(x, mean, sigma); } @@ -84,13 +83,12 @@ Int_t RooGaussian::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars double RooGaussian::analyticalIntegral(Int_t code, const char* rangeName) const { - using namespace RooFit::Detail::AnalyticalIntegrals; - if (code == 2) { - // Integration over mean if the code was "2" - return gaussianIntegral(mean.min(rangeName), mean.max(rangeName), x, sigma); - } - // Integration over x otherwise - return gaussianIntegral(x.min(rangeName), x.max(rangeName), mean, sigma); + using namespace RooFit::Detail::AnalyticalIntegrals; + + auto& constant = code == 1 ? mean : x; + auto& integrand = code == 1 ? x : mean; + + return gaussianIntegral(integrand.min(rangeName), integrand.max(rangeName), constant, sigma); } //////////////////////////////////////////////////////////////////////////////// @@ -136,11 +134,7 @@ void RooGaussian::generateEvent(Int_t code) void RooGaussian::translate(RooFit::Detail::CodeSquashContext &ctx) const { // Build a call to the stateless gaussian defined later. - std::string const& xName = ctx.getResult(&x.arg()); - std::string const& meanName = ctx.getResult(&mean.arg()); - std::string const& sigmaName = ctx.getResult(&sigma.arg()); - ctx.addResult(this, - "RooFit::Detail::EvaluateFuncs::gaussEvaluate(" + xName + ", " + meanName + ", " + sigmaName + ")"); + ctx.addResult(this, ctx.buildCall("RooFit::Detail::EvaluateFuncs::gaussianEvaluate", x, mean, sigma)); } //////////////////////////////////////////////////////////////////////////////// @@ -148,17 +142,9 @@ void RooGaussian::translate(RooFit::Detail::CodeSquashContext &ctx) const std::string RooGaussian::buildCallToAnalyticIntegral(Int_t code, const char *rangeName, RooFit::Detail::CodeSquashContext &ctx) const { - std::string const& sigmaName = ctx.getResult(&sigma.arg()); - if (code == 2) { - std::string const& meanMin = RooNumber::toString(mean.min(rangeName)); - std::string const& meanMax = RooNumber::toString(mean.max(rangeName)); - std::string const& xName = ctx.getResult(&x.arg()); - return "RooFit::Detail::AnalyticalIntegrals::gaussianIntegral(" + meanMin + ", " + meanMax + ", " + xName + ", " + - sigmaName + ")"; - } - std::string const& xMin = RooNumber::toString(x.min(rangeName)); - std::string const& xMax = RooNumber::toString(x.max(rangeName)); - std::string const& meanName = ctx.getResult(&mean.arg()); - return "RooFit::Detail::AnalyticalIntegrals::gaussianIntegral(" + xMin + ", " + xMax + ", " + meanName + ", " + - sigmaName + ")"; + auto& constant = code == 1 ? mean : x; + auto& integrand = code == 1 ? x : mean; + + return ctx.buildCall("RooFit::Detail::AnalyticalIntegrals::gaussianIntegral", + integrand.min(rangeName), integrand.max(rangeName), constant, sigma); } diff --git a/roofit/roofitcore/inc/RooFit/Detail/AnalyticalIntegrals.h b/roofit/roofitcore/inc/RooFit/Detail/AnalyticalIntegrals.h index c205dc9ee3015d4027d2d924133273f54cc04a19..e2f4ea613c5bf2a5b600aed0ba835a9591886758 100644 --- a/roofit/roofitcore/inc/RooFit/Detail/AnalyticalIntegrals.h +++ b/roofit/roofitcore/inc/RooFit/Detail/AnalyticalIntegrals.h @@ -16,6 +16,8 @@ #include <TMath.h> +#include <cmath> + namespace RooFit { namespace Detail { @@ -63,6 +65,15 @@ inline double gaussianIntegral(double xMin, double xMax, double mean, double sig return resultScale * 0.5 * cond; } +inline double exponentialIntegral(double xMin, double xMax, double constant) +{ + if (constant == 0.0) { + return xMax - xMin; + } + + return (std::exp(constant * xMax) - std::exp(constant * xMin)) / constant; +} + } // namespace AnalyticalIntegrals } // namespace Detail diff --git a/roofit/roofitcore/inc/RooFit/Detail/CodeSquashContext.h b/roofit/roofitcore/inc/RooFit/Detail/CodeSquashContext.h index bb0e55dbdc091b1a9c200f67d4c2b300c6da3644..36dff522948ebe918afb7a8a02ee2fd92325759a 100644 --- a/roofit/roofitcore/inc/RooFit/Detail/CodeSquashContext.h +++ b/roofit/roofitcore/inc/RooFit/Detail/CodeSquashContext.h @@ -14,11 +14,16 @@ #ifndef RooFit_Detail_CodeSquashContext_h #define RooFit_Detail_CodeSquashContext_h -#include "RooAbsArg.h" +#include <RooAbsArg.h> +#include <RooNumber.h> +#include <sstream> #include <string> #include <unordered_map> +template <class T> +class RooTemplateProxy; + namespace RooFit { namespace Detail { @@ -34,7 +39,13 @@ public: /// @brief Gets the result for the given node using the node name. /// @param key The node to get the result string for. /// @return String representing the result of this node. - inline std::string const &getResult(RooAbsArg const *key) const { return _nodeNames.at(key->namePtr()); } + inline std::string const &getResult(RooAbsArg const &key) const { return _nodeNames.at(key.namePtr()); } + + template <class T> + std::string const &getResult(RooTemplateProxy<T> const &key) const + { + return getResult(key.arg()); + } /// @brief Checks if the current node has a result string already assigned. /// @param key The node to get the result string for. @@ -69,7 +80,40 @@ public: /// @param in String to add to the squashed code. inline void addToCodeBody(std::string const &in) { _code += in; } + /// @brief Build the code to call the function with name `funcname`, passing some arguments. + /// The arguments can either be doubles or some RooFit arguments whose + /// results will be looked up in the context. + template <typename... Args_t> + std::string buildCall(std::string const &funcname, Args_t const &...args) const + { + std::stringstream ss; + ss << funcname << "(" << buildArgs(args...) << ")" << std::endl; + return ss.str(); + } + private: + std::string buildArg(double x) const { return RooNumber::toString(x); } + + template <class T> + std::string buildArg(T const &arg) const + { + return getResult(arg); + } + + std::string buildArgs() const { return ""; } + + template <class Arg_t> + std::string buildArgs(Arg_t const &arg) const + { + return buildArg(arg); + } + + template <typename Arg_t, typename... Args_t> + std::string buildArgs(Arg_t const &arg, Args_t const &...args) const + { + return buildArg(arg) + ", " + buildArgs(args...); + } + /// @brief Map of node names to their result strings. std::unordered_map<const TNamed *, std::string> _nodeNames; /// @brief Block of code that is placed before the rest of the function body. diff --git a/roofit/roofitcore/inc/RooFit/Detail/EvaluateFuncs.h b/roofit/roofitcore/inc/RooFit/Detail/EvaluateFuncs.h index 3b8d65e43ae13a45ecdc045647f3dc038da6e4e0..5f5848e76292a0802642d1912baf326ce769f9fb 100644 --- a/roofit/roofitcore/inc/RooFit/Detail/EvaluateFuncs.h +++ b/roofit/roofitcore/inc/RooFit/Detail/EvaluateFuncs.h @@ -23,7 +23,7 @@ namespace Detail { namespace EvaluateFuncs { /// @brief Function to evaluate an un-normalized RooGaussian. -inline double gaussEvaluate(double x, double mean, double sigma) +inline double gaussianEvaluate(double x, double mean, double sigma) { const double arg = x - mean; const double sig = sigma; diff --git a/roofit/roofitcore/inc/RooFuncWrapper.h b/roofit/roofitcore/inc/RooFuncWrapper.h index 9a5e9c8720d8697f24396b605d25c83401c3210a..13d8fe76ca78674b98d2e790c21af9ac5188d032 100644 --- a/roofit/roofitcore/inc/RooFuncWrapper.h +++ b/roofit/roofitcore/inc/RooFuncWrapper.h @@ -45,8 +45,8 @@ protected: double evaluate() const override; private: - std::string buildCode(RooAbsReal const &head, RooArgSet const & /* paramSet */, RooArgSet const &obsSet, - const RooAbsData *data); + std::string + buildCode(RooAbsReal const &head, RooArgSet const & /* paramSet */, RooArgSet const &obsSet, const RooAbsData *data); void updateGradientVarBuffer() const; diff --git a/roofit/roofitcore/src/RooAddition.cxx b/roofit/roofitcore/src/RooAddition.cxx index 8fd0d1ea1ac4768997a61b98ff1040891efa3aac..e4e2ced668004d39d6fe7ad59668d57904bbd45c 100644 --- a/roofit/roofitcore/src/RooAddition.cxx +++ b/roofit/roofitcore/src/RooAddition.cxx @@ -189,7 +189,7 @@ void RooAddition::translate(RooFit::Detail::CodeSquashContext &ctx) const std::string decl = "double " + varName + "[" + std::to_string(eleSize) + "]{"; int idx = 0; for (RooAbsArg *it : _set) { - decl += ctx.getResult(it) + ","; + decl += ctx.getResult(*it) + ","; ctx.addResult(it, varName + "[" + std::to_string(idx) + "]"); idx++; } @@ -208,7 +208,7 @@ void RooAddition::translate(RooFit::Detail::CodeSquashContext &ctx) const result = "("; for (RooAbsArg *it : _set) { - result += ctx.getResult(it) + '+'; + result += ctx.getResult(*it) + '+'; } result.back() = ')'; ctx.addResult(this, result); diff --git a/roofit/roofitcore/src/RooFuncWrapper.cxx b/roofit/roofitcore/src/RooFuncWrapper.cxx index 6aaa43045fed33f73cf16351e97a8b2632b5835c..2426cb5971547c661de99a9de3aa36ec77caebd7 100644 --- a/roofit/roofitcore/src/RooFuncWrapper.cxx +++ b/roofit/roofitcore/src/RooFuncWrapper.cxx @@ -43,11 +43,11 @@ RooFuncWrapper::RooFuncWrapper(const char *name, const char *title, RooAbsReal c std::unique_ptr<RooAbsReal> pdf{RooFit::Detail::compileForNormSet(obj, normSet)}; // Get the parameters. RooArgSet paramSet; - pdf->getParameters(data ? data->get() : nullptr, paramSet); + obj.getParameters(data ? data->get() : nullptr, paramSet); // Get the observable if we have a valid dataset. RooArgSet obsSet; if (data) - pdf->getObservables(data->get(), obsSet); + obj.getObservables(data->get(), obsSet); // Load the parameters and observables. loadParamsAndObs(name, paramSet, obsSet, data); @@ -58,7 +58,7 @@ RooFuncWrapper::RooFuncWrapper(const char *name, const char *title, RooAbsReal c declareAndDiffFunction(name, func); } -RooFuncWrapper::RooFuncWrapper(const RooFuncWrapper &other, const char *name /*=nullptr*/) +RooFuncWrapper::RooFuncWrapper(const RooFuncWrapper &other, const char *name) : RooAbsReal(other, name), _params("!params", this, other._params), _func(other._func), @@ -177,8 +177,8 @@ void RooFuncWrapper::gradient(const double *x, double *g) const _grad(const_cast<double *>(x), _observables.data(), g); } -std::string RooFuncWrapper::buildCode(RooAbsReal const &head, RooArgSet const & /* paramSet */, - RooArgSet const &obsSet, const RooAbsData *data) +std::string RooFuncWrapper::buildCode(RooAbsReal const &head, RooArgSet const & /* paramSet */, RooArgSet const &obsSet, + const RooAbsData *data) { RooFit::Detail::CodeSquashContext ctx; @@ -224,5 +224,5 @@ std::string RooFuncWrapper::buildCode(RooAbsReal const &head, RooArgSet const & curr->translate(ctx); } - return ctx.assembleCode(ctx.getResult(&head)); + return ctx.assembleCode(ctx.getResult(head)); } diff --git a/roofit/roofitcore/src/RooNormalizedPdf.cxx b/roofit/roofitcore/src/RooNormalizedPdf.cxx index cf668b00824625e774f285bac5578a24c2e6fe89..cc29d8877764b365592993596ae81017f5fc1ac7 100644 --- a/roofit/roofitcore/src/RooNormalizedPdf.cxx +++ b/roofit/roofitcore/src/RooNormalizedPdf.cxx @@ -20,7 +20,7 @@ */ void RooNormalizedPdf::computeBatch(cudaStream_t * /*stream*/, double *output, size_t nEvents, - RooFit::Detail::DataMap const& dataMap) const + RooFit::Detail::DataMap const &dataMap) const { auto nums = dataMap.at(_pdf); auto integralSpan = dataMap.at(_normIntegral); @@ -40,5 +40,5 @@ void RooNormalizedPdf::computeBatch(cudaStream_t * /*stream*/, double *output, s void RooNormalizedPdf::translate(RooFit::Detail::CodeSquashContext &ctx) const { // For now just return function/normalization integral. - ctx.addResult(this, ctx.getResult(&_pdf.arg()) + "/" + ctx.getResult(&_normIntegral.arg())); + ctx.addResult(this, ctx.getResult(_pdf) + "/" + ctx.getResult(_normIntegral)); } diff --git a/roofit/roofitcore/src/RooProduct.cxx b/roofit/roofitcore/src/RooProduct.cxx index e431cf4a5db1bfbc635195c2d6208c6fe418dfcc..30c0a23a37f3000788dc8554621aaa2d55d25cf1 100644 --- a/roofit/roofitcore/src/RooProduct.cxx +++ b/roofit/roofitcore/src/RooProduct.cxx @@ -500,8 +500,8 @@ void RooProduct::translate(RooFit::Detail::CodeSquashContext &ctx) const std::string result; // Build a (node1 * node2 * node3 * ...) like expression. result = '('; - for (const auto item : _compRSet) { - result += ctx.getResult(item) + "*"; + for (RooAbsArg* item : _compRSet) { + result += ctx.getResult(*item) + "*"; } result.back() = ')'; ctx.addResult(this, result); diff --git a/roofit/roofitcore/test/testRooFuncWrapper.cxx b/roofit/roofitcore/test/testRooFuncWrapper.cxx index d1e651c782bd63ea4aa34da102b66117b5b00035..5411446f127210df250c9f4ed356d85dddfbe8ef 100644 --- a/roofit/roofitcore/test/testRooFuncWrapper.cxx +++ b/roofit/roofitcore/test/testRooFuncWrapper.cxx @@ -16,12 +16,11 @@ #include <RooConstVar.h> #include <RooDataSet.h> #include <RooFitResult.h> +#include <RooExponential.h> #include <RooFuncWrapper.h> #include <RooGaussian.h> -#include <RooHelpers.h> #include <RooMinimizer.h> #include <RooProduct.h> -#include <RooRealIntegral.h> #include <RooRealVar.h> #include <TROOT.h> @@ -38,6 +37,12 @@ namespace { double getNumDerivative(const RooAbsReal &pdf, RooRealVar &var, const RooArgSet &normSet, double eps = 1e-8) { double orig = var.getVal(); + if (!var.inRange(orig + eps, nullptr)) { + throw std::runtime_error("getNumDerivative(): positive variation outside of range!"); + } + if (!var.inRange(orig - eps, nullptr)) { + throw std::runtime_error("getNumDerivative(): negative variation outside of range!"); + } var.setVal(orig + eps); double plus = pdf.getVal(normSet); var.setVal(orig - eps); @@ -193,8 +198,6 @@ TEST(RooFuncWrapper, NllWithObservables) TEST(RooFuncWrapper, GaussianNormalized) { - using namespace RooFit; - RooRealVar x("x", "x", 0, -10, std::numeric_limits<double>::infinity()); RooRealVar mu("mu", "mu", 0, -10, 10); @@ -229,3 +232,30 @@ TEST(RooFuncWrapper, GaussianNormalized) EXPECT_NEAR(getNumDerivative(gauss, static_cast<RooRealVar &>(*paramsGauss[i]), normSet), dMyGauss[i], 1e-8); } } + +TEST(RooFuncWrapper, Exponential) +{ + RooRealVar x("x", "x", 1.0, 0, 10); + RooRealVar c("c", "c", 0.1, 0, 10); + + RooExponential expo("expo", "expo", x, c); + + RooArgSet normSet{x}; + + RooFuncWrapper expoFunc("expo", "expo", expo, normSet); + + RooArgSet params; + expo.getParameters(nullptr, params); + + EXPECT_NEAR(expo.getVal(normSet), expoFunc.getVal(), 1e-8); + + // Get AD based derivative + double dExpo[2] = {}; + expoFunc.getGradient(dExpo); + + // Check if derivatives are equal + for (std::size_t i = 0; i < params.size(); ++i) { + EXPECT_NEAR(getNumDerivative(expo, static_cast<RooRealVar &>(*params[i]), normSet), dExpo[i], 1e-8) + << params[i]->GetName(); + } +}