diff --git a/interpreter/cling/tools/plugins/clad/CMakeLists.txt b/interpreter/cling/tools/plugins/clad/CMakeLists.txt index 79e1cd41d23d1..f0d3d843397ba 100644 --- a/interpreter/cling/tools/plugins/clad/CMakeLists.txt +++ b/interpreter/cling/tools/plugins/clad/CMakeLists.txt @@ -74,7 +74,7 @@ if (DEFINED CLAD_SOURCE_DIR) list(APPEND _clad_extra_settings SOURCE_DIR ${CLAD_SOURCE_DIR}) else() list(APPEND _clad_extra_settings GIT_REPOSITORY https://github.com/vgvassilev/clad.git) - list(APPEND _clad_extra_settings GIT_TAG v2.2) + list(APPEND _clad_extra_settings GIT_TAG master) endif() ## list(APPEND _clad_patches_list "patch1.patch" "patch2.patch") diff --git a/roofit/codegen/src/CodegenImpl.cxx b/roofit/codegen/src/CodegenImpl.cxx index 4fd12803b13b6..e0c3cefb42b05 100644 --- a/roofit/codegen/src/CodegenImpl.cxx +++ b/roofit/codegen/src/CodegenImpl.cxx @@ -106,7 +106,7 @@ void rooHistTranslateImpl(RooAbsArg const &arg, CodegenContext &ctx, int intOrde } std::string realSumPdfTranslateImpl(CodegenContext &ctx, RooAbsArg const &arg, RooArgList const &funcList, - RooArgList const &coefList, bool normalize) + RooArgList const &coefList, bool normalize, bool forceScopeIndependent) { bool noLastCoeff = funcList.size() != coefList.size(); @@ -116,7 +116,12 @@ std::string realSumPdfTranslateImpl(CodegenContext &ctx, RooAbsArg const &arg, R std::string sum = ctx.getTmpVarName(); std::string coeffSum = ctx.getTmpVarName(); - ctx.addToCodeBody(&arg, "double " + sum + " = 0;\ndouble " + coeffSum + "= 0;\n"); + std::string code1 = "double " + sum + " = 0;\ndouble " + coeffSum + "= 0;\n"; + + if (forceScopeIndependent) + ctx.addToCodeBody(code1, true); + else + ctx.addToCodeBody(&arg, code1); std::string iterator = "i_" + ctx.getTmpVarName(); std::string subscriptExpr = "[" + iterator + "]"; @@ -131,7 +136,10 @@ std::string realSumPdfTranslateImpl(CodegenContext &ctx, RooAbsArg const &arg, R } else if (normalize) { code += sum + " /= " + coeffSum + ";\n"; } - ctx.addToCodeBody(&arg, code); + if (forceScopeIndependent) + ctx.addToCodeBody(code, true); + else + ctx.addToCodeBody(&arg, code); return sum; } @@ -243,7 +251,7 @@ void codegenImpl(RooAbsArg &arg, CodegenContext &ctx) void codegenImpl(RooAddPdf &arg, CodegenContext &ctx) { - ctx.addResult(&arg, realSumPdfTranslateImpl(ctx, arg, arg.pdfList(), arg.coefList(), true)); + ctx.addResult(&arg, realSumPdfTranslateImpl(ctx, arg, arg.pdfList(), arg.coefList(), true, false)); } void codegenImpl(RooMultiVarGaussian &arg, CodegenContext &ctx) @@ -264,30 +272,25 @@ void codegenImpl(RooMultiPdf &arg, CodegenContext &ctx) // indices MathFunc call becomes more efficient. if (numPdfs > 2) { ctx.addResult(&arg, ctx.buildCall(mathFunc("multipdf"), arg.indexCategory(), arg.getPdfList())); + return; + } + // Ternary nested expression + std::string indexExpr = ctx.getResult(arg.indexCategory()); - std::cout << "MathFunc call used\n"; - - } else { - - // Ternary nested expression - std::string indexExpr = ctx.getResult(arg.indexCategory()); - - // int numPdfs = arg.getNumPdfs(); - std::string expr; + // int numPdfs = arg.getNumPdfs(); + std::string expr; - for (int i = 0; i < numPdfs; ++i) { - RooAbsPdf *pdf = arg.getPdf(i); - std::string pdfExpr = ctx.getResult(*pdf); + for (int i = 0; i < numPdfs; ++i) { + RooAbsPdf *pdf = arg.getPdf(i); + std::string pdfExpr = ctx.getResult(*pdf); - expr += "(" + indexExpr + " == " + std::to_string(i) + " ? (" + pdfExpr + ") : "; - } + expr += "(" + indexExpr + " == " + std::to_string(i) + " ? (" + pdfExpr + ") : "; + } - expr += "0.0"; - expr += std::string(numPdfs, ')'); // Close all ternary operators + expr += "0.0"; + expr += std::string(numPdfs, ')'); // Close all ternary operators - ctx.addResult(&arg, expr); - std::cout << "Ternary expression call used \n"; - } + ctx.addResult(&arg, expr); } // RooCategory index added. @@ -308,6 +311,7 @@ void codegenImpl(RooAddition &arg, CodegenContext &ctx) { if (arg.list().empty()) { ctx.addResult(&arg, "0.0"); + return; } std::string result; if (arg.list().size() > 1) @@ -558,7 +562,6 @@ void codegenImpl(RooFit::Detail::RooNLLVarNew &arg, CodegenContext &ctx) std::string weightSumName = RooFit::Detail::makeValidVarName(arg.GetName()) + "WeightSum"; std::string resName = RooFit::Detail::makeValidVarName(arg.GetName()) + "Result"; - ctx.addResult(&arg, resName); ctx.addToGlobalScope("double " + weightSumName + " = 0.0;\n"); ctx.addToGlobalScope("double " + resName + " = 0.0;\n"); @@ -585,6 +588,8 @@ void codegenImpl(RooFit::Detail::RooNLLVarNew &arg, CodegenContext &ctx) std::string expected = ctx.getResult(*arg.expectedEvents()); ctx.addToCodeBody(resName + " += " + expected + " - " + weightSumName + " * std::log(" + expected + ");\n"); } + + ctx.addResult(&arg, resName); } void codegenImpl(RooFit::Detail::RooNormalizedPdf &arg, CodegenContext &ctx) @@ -642,7 +647,14 @@ void codegenImpl(RooPolynomial &arg, CodegenContext &ctx) void codegenImpl(RooProduct &arg, CodegenContext &ctx) { - ctx.addResult(&arg, ctx.buildCall(mathFunc("product"), arg.realComponents(), arg.realComponents().size())); + std::stringstream ss; + std::size_t n = arg.realComponents().size(); + for (std::size_t i = 0; i < n; ++i) { + ss << ctx.getResult(arg.realComponents()[i]); + if (i != n - 1) + ss << " * " << std::endl; + } + ctx.addResult(&arg, ss.str()); } void codegenImpl(RooRatio &arg, CodegenContext &ctx) @@ -698,17 +710,17 @@ void codegenImpl(RooRealIntegral &arg, CodegenContext &ctx) auto &intVar = static_cast(*arg.numIntRealVars()[0]); std::string obsName = ctx.getTmpVarName(); - std::string oldIntVarResult = ctx.getResult(intVar); - ctx.addResult(&intVar, "obs[0]"); + auto oldVecObsInfo = ctx._vecObsIndices[intVar.namePtr()]; + ctx.addVecObs(intVar.GetName(), 0); std::string funcName = ctx.buildFunction(arg.integrand(), {}); + ctx._vecObsIndices[intVar.namePtr()] = oldVecObsInfo; std::stringstream ss; ss << "double " << obsName << "[1];\n"; std::string resName = RooFit::Detail::makeValidVarName(arg.GetName()) + "Result"; - ctx.addResult(&arg, resName); ctx.addToGlobalScope("double " + resName + " = 0.0;\n"); // TODO: once Clad has support for higher-order functions (follow also the @@ -730,24 +742,21 @@ void codegenImpl(RooRealIntegral &arg, CodegenContext &ctx) ctx.addToGlobalScope(ss.str()); - ctx.addResult(&intVar, oldIntVarResult); + ctx.addResult(&arg, resName); } void codegenImpl(RooRealSumFunc &arg, CodegenContext &ctx) { - ctx.addResult(&arg, realSumPdfTranslateImpl(ctx, arg, arg.funcList(), arg.coefList(), false)); + ctx.addResult(&arg, realSumPdfTranslateImpl(ctx, arg, arg.funcList(), arg.coefList(), false, false)); } void codegenImpl(RooRealSumPdf &arg, CodegenContext &ctx) { - ctx.addResult(&arg, realSumPdfTranslateImpl(ctx, arg, arg.funcList(), arg.coefList(), false)); + ctx.addResult(&arg, realSumPdfTranslateImpl(ctx, arg, arg.funcList(), arg.coefList(), false, false)); } void codegenImpl(RooRealVar &arg, CodegenContext &ctx) { - if (!arg.isConstant()) { - ctx.addResult(&arg, arg.GetName()); - } ctx.addResult(&arg, doubleToString(arg.getVal())); } @@ -988,7 +997,7 @@ std::string codegenIntegralImpl(RooPolynomial &arg, int, const char *rangeName, std::string codegenIntegralImpl(RooRealSumPdf &arg, int code, const char *rangeName, CodegenContext &ctx) { // Re-use translate, since integration is linear. - return realSumPdfTranslateImpl(ctx, arg, arg.funcIntListFromCache(code, rangeName), arg.coefList(), false); + return realSumPdfTranslateImpl(ctx, arg, arg.funcIntListFromCache(code, rangeName), arg.coefList(), false, true); } std::string codegenIntegralImpl(RooUniform &arg, int code, const char *rangeName, CodegenContext &) diff --git a/roofit/roofitcore/inc/RooFit/CodegenContext.h b/roofit/roofitcore/inc/RooFit/CodegenContext.h index f61b3f08450bb..0d7b81edd4b86 100644 --- a/roofit/roofitcore/inc/RooFit/CodegenContext.h +++ b/roofit/roofitcore/inc/RooFit/CodegenContext.h @@ -30,8 +30,7 @@ template class RooTemplateProxy; -namespace RooFit { -namespace Experimental { +namespace RooFit::Experimental { template struct Prio { @@ -46,18 +45,18 @@ using PrioLowest = Prio<10>; class CodegenContext { public: void addResult(RooAbsArg const *key, std::string const &value); - void addResult(const char *key, std::string const &value); - std::string const &getResult(RooAbsArg const &arg); + std::string getResult(RooAbsArg const &arg); template - std::string const &getResult(RooTemplateProxy const &key) + std::string getResult(RooTemplateProxy const &key) { return getResult(key.arg()); } void addToGlobalScope(std::string const &str); void addVecObs(const char *key, int idx); + void addParam(const RooAbsArg *key, int idx); int observableIndexOf(const RooAbsArg &arg) const; void addToCodeBody(RooAbsArg const *klass, std::string const &in); @@ -123,6 +122,13 @@ class CodegenContext { }; ScopeRAII OutputScopeRangeComment(RooAbsArg const *arg) { return {arg, *this}; } + /// @brief Map of node names to their result strings. + std::unordered_map _nodeNames; + std::size_t _nWksp = 0; + std::unordered_map _paramIndices; + /// @brief A map to keep track of the observable indices if they are non scalar. + std::unordered_map _vecObsIndices; + private: void pushScope(); void popScope(); @@ -133,8 +139,6 @@ class CodegenContext { void endLoop(LoopScope const &scope); - void addResult(TNamed const *key, std::string const &value); - template {}, bool>::type = true> std::string buildArg(T x) { @@ -179,10 +183,6 @@ class CodegenContext { template std::string typeName() const; - /// @brief Map of node names to their result strings. - std::unordered_map _nodeNames; - /// @brief A map to keep track of the observable indices if they are non scalar. - std::unordered_map _vecObsIndices; /// @brief Indicate whether a node depends on the dataset. std::unordered_set _dependsOnData; /// @brief The code layered by lexical scopes used as a stack. @@ -191,8 +191,6 @@ class CodegenContext { unsigned _indent = 0; /// @brief Index to get unique names for temporary variables. mutable int _tmpVarIdx = 0; - /// @brief A map to keep track of list names as assigned by addResult. - std::unordered_map::Value_t, std::string> _listNames; std::vector _xlArr; std::vector _collectedFunctions; std::string _collectedCode; @@ -231,7 +229,6 @@ void declareDispatcherCode(std::string const &funcName); void codegen(RooAbsArg &arg, CodegenContext &ctx); -} // namespace Experimental -} // namespace RooFit +} // namespace RooFit::Experimental #endif diff --git a/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h b/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h index e9cef8def1f61..64d8245e1d9a4 100644 --- a/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h +++ b/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h @@ -820,7 +820,23 @@ double stepFunctionIntegral(double xmin, double xmax, std::size_t nBins, DoubleA } // namespace RooFit::Detail::MathFuncs +inline void fillFromWorkspace(double *out, std::size_t n, double const *wksp, double const *idx) +{ + for (std::size_t i = 0; i < n; ++i) { + out[i] += wksp[static_cast(idx[i])]; + } +} + namespace clad::custom_derivatives { + +inline void fillFromWorkspace_pullback(double *, std::size_t n, double const *, double const *idx, double *d_out, + std::size_t *, double *d_wksp, double *) +{ + for (std::size_t i = 0; i < n; ++i) { + d_wksp[static_cast(idx[i])] += d_out[i]; + } +} + namespace RooFit::Detail::MathFuncs { // Clad can't generate the pullback for binNumber because of the @@ -834,6 +850,7 @@ void binNumber_pullback(Types...) } } // namespace RooFit::Detail::MathFuncs + } // namespace clad::custom_derivatives #endif diff --git a/roofit/roofitcore/src/RooEvaluatorWrapper.cxx b/roofit/roofitcore/src/RooEvaluatorWrapper.cxx index 9f4d3be01b532..27fd38e81c6cf 100644 --- a/roofit/roofitcore/src/RooEvaluatorWrapper.cxx +++ b/roofit/roofitcore/src/RooEvaluatorWrapper.cxx @@ -250,13 +250,12 @@ RooFuncWrapper::RooFuncWrapper(RooAbsReal &obj, const RooAbsData *data, RooSimul // First update the result variable of params in the compute graph to in[]. int idx = 0; for (RooAbsArg *param : _params) { - ctx.addResult(param, "params[" + std::to_string(idx) + "]"); + ctx.addParam(param, idx); idx++; } for (auto const &item : _obsInfos) { const char *obsName = item.first->GetName(); - ctx.addResult(obsName, "obs"); ctx.addVecObs(obsName, item.second); } diff --git a/roofit/roofitcore/src/RooFit/CodegenContext.cxx b/roofit/roofitcore/src/RooFit/CodegenContext.cxx index 8a1f10e7c958e..76d0bcc132923 100644 --- a/roofit/roofitcore/src/RooFit/CodegenContext.cxx +++ b/roofit/roofitcore/src/RooFit/CodegenContext.cxx @@ -25,39 +25,15 @@ #include #include -namespace { - -bool startsWith(std::string_view str, std::string_view prefix) -{ - return str.size() >= prefix.size() && 0 == str.compare(0, prefix.size(), prefix); -} - -} // namespace - namespace RooFit { namespace Experimental { -/// @brief Adds (or overwrites) the string representing the result of a node. -/// @param key The name of the node to add the result for. -/// @param value The new name to assign/overwrite. -void CodegenContext::addResult(const char *key, std::string const &value) -{ - const TNamed *namePtr = RooNameReg::known(key); - if (namePtr) - addResult(namePtr, value); -} - -void CodegenContext::addResult(TNamed const *key, std::string const &value) -{ - _nodeNames[key] = value; -} - /// @brief Gets the result for the given node using the node name. This node also performs the necessary /// code generation through recursive calls to 'translate'. A call to this function modifies the already /// existing code body. /// @param key The node to get the result string for. /// @return String representing the result of this node. -std::string const &CodegenContext::getResult(RooAbsArg const &arg) +std::string CodegenContext::getResult(RooAbsArg const &arg) { // If the result has already been recorded, just return the result. // It is usually the responsibility of each translate function to assign @@ -65,23 +41,19 @@ std::string const &CodegenContext::getResult(RooAbsArg const &arg) // for a particular node, it means the node has already been 'translate'd and we // dont need to visit it again. auto found = _nodeNames.find(arg.namePtr()); - if (found != _nodeNames.end()) - return found->second; - - // The result for vector observables should already be in the map if you - // opened the loop scope. This is just to check if we did not request the - // result of a vector-valued observable outside of the scope of a loop. - auto foundVecObs = _vecObsIndices.find(arg.namePtr()); - if (foundVecObs != _vecObsIndices.end()) { - throw std::runtime_error("You requested the result of a vector observable outside a loop scope for it!"); - } + std::size_t idx = 0; + if (found != _nodeNames.end()) { + idx = found->second; + } else { - auto RAII(OutputScopeRangeComment(&arg)); + auto RAII(OutputScopeRangeComment(&arg)); - // Now, recursively call translate into the current argument to load the correct result. - codegen(const_cast(arg), *this); + // Now, recursively call translate into the current argument to load the correct result. + codegen(const_cast(arg), *this); - return _nodeNames.at(arg.namePtr()); + idx = _nodeNames.at(arg.namePtr()); + } + return "wksp[" + std::to_string(idx) + "]"; } /// @brief Adds the given string to the string block that will be emitted at the top of the squashed function. Useful @@ -105,6 +77,11 @@ void CodegenContext::addVecObs(const char *key, int idx) _vecObsIndices[namePtr] = idx; } +void CodegenContext::addParam(RooAbsArg const *arg, int idx) +{ + _paramIndices[arg] = idx; +} + int CodegenContext::observableIndexOf(RooAbsArg const &arg) const { auto it = _vecObsIndices.find(arg.namePtr()); @@ -175,7 +152,6 @@ std::unique_ptr CodegenContext::beginLoop(RooAbsArg c continue; vars.push_back(it.first); - _nodeNames[it.first] = "obs[static_cast(obs[" + std::to_string(2 * it.second) + "]) + " + idx + "]"; if (firstObsIdx == -1) { firstObsIdx = it.second; } @@ -190,6 +166,18 @@ std::unique_ptr CodegenContext::beginLoop(RooAbsArg c addToCodeBody(in, "for(int " + idx + " = 0; " + idx + " < obs[" + std::to_string(2 * firstObsIdx + 1) + "]; " + idx + "++) {\n"); + // set the results of the vector observables + for (auto const &it : _vecObsIndices) { + if (!in->dependsOn(it.first)) + continue; + + auto savedName = "wksp[" + std::to_string(_nWksp) + "]"; + std::string outVarDecl; + outVarDecl = savedName + " = obs[static_cast(obs[" + std::to_string(2 * it.second) + "]) + " + idx + "];\n"; + addToCodeBody(outVarDecl); + _nodeNames[it.first] = _nWksp++; + } + return std::make_unique(*this, std::move(vars)); } @@ -216,28 +204,8 @@ std::string CodegenContext::getTmpVarName() const /// @param valueToSave The actual string value to save as a temporary. void CodegenContext::addResult(RooAbsArg const *in, std::string const &valueToSave) { - // std::string savedName = RooFit::Detail::makeValidVarName(in->GetName()); - std::string savedName = getTmpVarName(); - - // Only save values if they contain operations or they are numerals. Otherwise, we can use them directly. - - // Check if string is numeric. - char *end; - std::strtod(valueToSave.c_str(), &end); - bool isNumeric = (*end == '\0'); - - const bool hasOperations = valueToSave.find_first_of(":-+/*") != std::string::npos; - - // If the name is not empty and this value is worth saving, save it to the correct scope. - // otherwise, just return the actual value itself - if (hasOperations || isNumeric) { - std::string outVarDecl = "const double " + savedName + " = " + valueToSave + ";\n"; - addToCodeBody(in, outVarDecl); - } else { - savedName = valueToSave; - } - - addResult(in->namePtr(), savedName); + addToCodeBody(in, "wksp[" + std::to_string(_nWksp) + "]" + " = " + valueToSave + ";\n"); + _nodeNames[in->namePtr()] = _nWksp++; } /// @brief Function to save a RooListProxy as an array in the squashed code. @@ -249,25 +217,24 @@ std::string CodegenContext::buildArg(RooAbsCollection const &in, std::string con return "nullptr"; } - auto it = _listNames.find(in.uniqueId().value()); - if (it != _listNames.end()) - return it->second; - std::string savedName = getTmpVarName(); bool canSaveOutside = true; + std::vector indices; + indices.reserve(in.size()); + std::stringstream declStrm; - declStrm << arrayType << " " << savedName << "[]{"; + declStrm << arrayType << " " << savedName << "[" << in.size() << "]{};\n"; for (const auto arg : in) { - declStrm << getResult(*arg) << ","; + getResult(*arg); // fill the result cache + indices.push_back(_nodeNames.at(arg->namePtr())); canSaveOutside = canSaveOutside && isScopeIndependent(arg); } - declStrm.seekp(-1, declStrm.cur); - declStrm << "};\n"; + + declStrm << "fillFromWorkspace(" << savedName << ", " << in.size() << ", wksp, " << buildArg(indices) << ");\n"; addToCodeBody(declStrm.str(), canSaveOutside); - _listNames.insert({in.uniqueId().value(), savedName}); return savedName; } @@ -334,12 +301,7 @@ CodegenContext::buildFunction(RooAbsArg const &arg, std::unordered_set::infinity();\n" - << funcBody << "\n}\n\n"; + << "constexpr double inf = std::numeric_limits::infinity();\n"; + if (ctx._nWksp > 0) { + bodyWithSigStrm << "double wksp [" << ctx._nWksp << "]{};\n"; + } + bodyWithSigStrm << funcBody << "\n}"; ctx._collectedFunctions.emplace_back(funcName); ctx._collectedCode += bodyWithSigStrm.str(); @@ -402,6 +367,20 @@ struct Caller_FUNC_NAME { void codegen(RooAbsArg &arg, CodegenContext &ctx) { + // parameters + auto foundParam = ctx._paramIndices.find(&arg); + if (foundParam != ctx._paramIndices.end()) { + ctx.addResult(&arg, "params[" + std::to_string(foundParam->second) + "]"); + return; + } + + // observables + auto foundObs = ctx._vecObsIndices.find(arg.namePtr()); + if (foundObs != ctx._vecObsIndices.end()) { + ctx.addResult(&arg, "obs[" + std::to_string(foundObs->second) + "]"); + return; + } + static bool codeDeclared = false; if (!codeDeclared) { declareDispatcherCode("codegenImpl");