Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooFuncWrapper.cxx
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 * Authors:
4 * Garima Singh, CERN 2022
5 *
6 * Copyright (c) 2022, CERN
7 *
8 * Redistribution and use in source and binary forms,
9 * with or without modification, are permitted according to the terms
10 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
11 */
12
13#include <RooFuncWrapper.h>
14
15#include <RooAbsData.h>
18#include <RooFit/Evaluator.h>
19#include <RooGlobalFunc.h>
20#include <RooHelpers.h>
21#include <RooMsgService.h>
22#include <RooRealVar.h>
23#include <RooSimultaneous.h>
24#include "RooEvaluatorWrapper.h"
25
26#include <TROOT.h>
27#include <TSystem.h>
28
29#include <fstream>
30#include <set>
31
32namespace {
33
34void replaceAll(std::string &str, const std::string &from, const std::string &to)
35{
36 if (from.empty())
37 return;
38 size_t start_pos = 0;
39 while ((start_pos = str.find(from, start_pos)) != std::string::npos) {
40 str.replace(start_pos, from.length(), to);
41 start_pos += to.length(); // In case 'to' contains 'from', like replacing 'x' with 'yx'
42 }
43}
44
45} // namespace
46
47namespace RooFit {
48
49namespace Experimental {
50
51RooFuncWrapper::RooFuncWrapper(const char *name, const char *title, RooAbsReal &obj, const RooAbsData *data,
52 RooSimultaneous const *simPdf, bool useEvaluator)
53 : RooAbsReal{name, title}, _params{"!params", "List of parameters", this}, _useEvaluator{useEvaluator}
54{
55 if (_useEvaluator) {
56 _absReal = std::make_unique<RooEvaluatorWrapper>(obj, const_cast<RooAbsData *>(data), false, "", simPdf, false);
57 }
58
59 std::string func;
60
61 // Get the parameters.
62 RooArgSet paramSet;
63 obj.getParameters(data ? data->get() : nullptr, paramSet);
64 RooArgSet floatingParamSet;
65 for (RooAbsArg *param : paramSet) {
66 if (!param->isConstant()) {
67 floatingParamSet.add(*param);
68 }
69 }
70
71 // Load the parameters and observables.
72 loadParamsAndData(&obj, floatingParamSet, data, simPdf);
73
74 func = buildCode(obj);
75
76 gInterpreter->Declare("#pragma cling optimize(2)");
77
78 // Declare the function and create its derivative.
80 _func = reinterpret_cast<Func>(gInterpreter->ProcessLine((_funcName + ";").c_str()));
81}
82
84 : RooAbsReal(other, name),
85 _params("!params", this, other._params),
86 _funcName(other._funcName),
87 _func(other._func),
88 _grad(other._grad),
89 _hasGradient(other._hasGradient),
90 _gradientVarBuffer(other._gradientVarBuffer),
91 _observables(other._observables)
92{
93}
94
95void RooFuncWrapper::loadParamsAndData(RooAbsArg const *head, RooArgSet const &paramSet, const RooAbsData *data,
96 RooSimultaneous const *simPdf)
97{
98 // Extract observables
99 std::stack<std::vector<double>> vectorBuffers; // for data loading
100 std::map<RooFit::Detail::DataKey, std::span<const double>> spans;
101
102 if (data) {
103 spans = RooFit::Detail::BatchModeDataHelpers::getDataSpans(*data, "", simPdf, true, false, vectorBuffers);
104 }
105
106 std::size_t idx = 0;
107 for (auto const &item : spans) {
108 std::size_t n = item.second.size();
109 _obsInfos.emplace(item.first, ObsInfo{idx, n});
110 _observables.reserve(_observables.size() + n);
111 for (std::size_t i = 0; i < n; ++i) {
112 _observables.push_back(item.second[i]);
113 }
114 idx += n;
115 }
116
117 // Extract parameters
118 for (auto *param : paramSet) {
119 if (!dynamic_cast<RooAbsReal *>(param)) {
120 std::stringstream errorMsg;
121 errorMsg << "In creation of function " << GetName()
122 << " wrapper: input param expected to be of type RooAbsReal.";
123 coutE(InputArguments) << errorMsg.str() << std::endl;
124 throw std::runtime_error(errorMsg.str().c_str());
125 }
126 if (spans.find(param) == spans.end()) {
127 _params.add(*param);
128 }
129 }
130 _gradientVarBuffer.resize(_params.size());
131
132 if (head) {
133 _nodeOutputSizes = RooFit::Detail::BatchModeDataHelpers::determineOutputSizes(
134 *head, [&spans](RooFit::Detail::DataKey key) -> int {
135 auto found = spans.find(key);
136 return found != spans.end() ? found->second.size() : -1;
137 });
138 }
139}
140
141std::string RooFuncWrapper::declareFunction(std::string const &funcBody)
142{
143 static int iFuncWrapper = 0;
144 auto funcName = "roo_func_wrapper_" + std::to_string(iFuncWrapper++);
145
146 // Declare the function
147 std::stringstream bodyWithSigStrm;
148 bodyWithSigStrm << "double " << funcName << "(double* params, double const* obs, double const* xlArr) {\n"
149 << funcBody << "\n}";
150 _collectedFunctions.emplace_back(funcName);
151 if (!gInterpreter->Declare(bodyWithSigStrm.str().c_str())) {
152 std::stringstream errorMsg;
153 errorMsg << "Function " << funcName << " could not be compiled. See above for details.";
154 oocoutE(nullptr, InputArguments) << errorMsg.str() << std::endl;
155 throw std::runtime_error(errorMsg.str().c_str());
156 }
157 return funcName;
158}
159
160void RooFuncWrapper::createGradient()
161{
162 std::string gradName = _funcName + "_grad_0";
163 std::string requestName = _funcName + "_req";
164
165 // Calculate gradient
166 gInterpreter->Declare("#include <Math/CladDerivator.h>\n");
167 // disable clang-format for making the following code unreadable.
168 // clang-format off
169 std::stringstream requestFuncStrm;
170 requestFuncStrm << "#pragma clad ON\n"
171 "void " << requestName << "() {\n"
172 " clad::gradient(" << _funcName << ", \"params\");\n"
173 "}\n"
174 "#pragma clad OFF";
175 // clang-format on
176 if (!gInterpreter->Declare(requestFuncStrm.str().c_str())) {
177 std::stringstream errorMsg;
178 errorMsg << "Function " << GetName() << " could not be differentiated. See above for details.";
179 oocoutE(nullptr, InputArguments) << errorMsg.str() << std::endl;
180 throw std::runtime_error(errorMsg.str().c_str());
181 }
182
183 _grad = reinterpret_cast<Grad>(gInterpreter->ProcessLine((gradName + ";").c_str()));
184 _hasGradient = true;
185}
186
187void RooFuncWrapper::gradient(double *out) const
188{
189 updateGradientVarBuffer();
190 std::fill(out, out + _params.size(), 0.0);
191
192 _grad(_gradientVarBuffer.data(), _observables.data(), _xlArr.data(), out);
193}
194
195void RooFuncWrapper::updateGradientVarBuffer() const
196{
197 std::transform(_params.begin(), _params.end(), _gradientVarBuffer.begin(),
198 [](RooAbsArg *obj) { return static_cast<RooAbsReal *>(obj)->getVal(); });
199}
200
201double RooFuncWrapper::evaluate() const
202{
203 if (_useEvaluator)
204 return _absReal->getVal();
205 updateGradientVarBuffer();
206
207 return _func(_gradientVarBuffer.data(), _observables.data(), _xlArr.data());
208}
209
210void RooFuncWrapper::gradient(const double *x, double *g) const
211{
212 std::fill(g, g + _params.size(), 0.0);
213
214 _grad(const_cast<double *>(x), _observables.data(), _xlArr.data(), g);
215}
216
217std::string RooFuncWrapper::buildCode(RooAbsReal const &head)
218{
219 RooFit::Detail::CodeSquashContext ctx(_nodeOutputSizes, _xlArr, *this);
220
221 // First update the result variable of params in the compute graph to in[<position>].
222 int idx = 0;
223 for (RooAbsArg *param : _params) {
224 ctx.addResult(param, "params[" + std::to_string(idx) + "]");
225 idx++;
226 }
227
228 for (auto const &item : _obsInfos) {
229 const char *name = item.first->GetName();
230 // If the observable is scalar, set name to the start idx. else, store
231 // the start idx and later set the the name to obs[start_idx + curr_idx],
232 // here curr_idx is defined by a loop producing parent node.
233 if (item.second.size == 1) {
234 ctx.addResult(name, "obs[" + std::to_string(item.second.idx) + "]");
235 } else {
236 ctx.addResult(name, "obs");
237 ctx.addVecObs(name, item.second.idx);
238 }
239 }
240
241 return ctx.assembleCode(ctx.getResult(head));
242}
243
244/// @brief Dumps a macro "filename.C" that can be used to test and debug the generated code and gradient.
245void RooFuncWrapper::writeDebugMacro(std::string const &filename) const
246{
247 std::stringstream allCode;
248 std::set<std::string> seenFunctions;
249
250 // Remove duplicated declared functions
251 for (std::string const &name : _collectedFunctions) {
252 if (seenFunctions.count(name) > 0) {
253 continue;
254 }
255 seenFunctions.insert(name);
256 std::unique_ptr<TInterpreterValue> v = gInterpreter->MakeInterpreterValue();
257 gInterpreter->Evaluate(name.c_str(), *v);
258 std::string s = v->ToString();
259 for (int i = 0; i < 2; ++i) {
260 s = s.erase(0, s.find("\n") + 1);
261 }
262 allCode << s << std::endl;
263 }
264
265 std::ofstream outFile;
266 outFile.open(filename + ".C");
267 outFile << R"(//auto-generated test macro
268#include <RooFit/Detail/MathFuncs.h>
269#include <Math/CladDerivator.h>
270
271#pragma cling optimize(2)
272)" << allCode.str()
273 << R"(
274#pragma clad ON
275void gradient_request() {
276 clad::gradient()"
277 << _funcName << R"(, "params");
278}
279#pragma clad OFF
280)";
281
282 updateGradientVarBuffer();
283
284 auto writeVector = [&](std::string const &name, std::span<const double> vec) {
285 std::stringstream decl;
286 decl << "std::vector<double> " << name << " = {";
287 for (std::size_t i = 0; i < vec.size(); ++i) {
288 if (i % 10 == 0)
289 decl << "\n ";
290 decl << vec[i];
291 if (i < vec.size() - 1)
292 decl << ", ";
293 }
294 decl << "\n};\n";
295
296 std::string declStr = decl.str();
297
298 replaceAll(declStr, "inf", "std::numeric_limits<double>::infinity()");
299 replaceAll(declStr, "nan", "NAN");
300
301 outFile << declStr;
302 };
303
304 outFile << "// clang-format off\n" << std::endl;
305 writeVector("parametersVec", _gradientVarBuffer);
306 outFile << std::endl;
307 writeVector("observablesVec", _observables);
308 outFile << std::endl;
309 writeVector("auxConstantsVec", _xlArr);
310 outFile << std::endl;
311 outFile << "// clang-format on\n" << std::endl;
312
313 outFile << R"(
314// To run as a ROOT macro
315void )" << filename
316 << R"(()
317{
318 std::vector<double> gradientVec(parametersVec.size());
319
320 auto func = [&](std::span<double> params) {
321 return )"
322 << _funcName << R"((params.data(), observablesVec.data(), auxConstantsVec.data());
323 };
324 auto grad = [&](std::span<double> params, std::span<double> out) {
325 return )"
326 << _funcName << R"(_grad_0(parametersVec.data(), observablesVec.data(), auxConstantsVec.data(),
327 out.data());
328 };
329
330 grad(parametersVec, gradientVec);
331
332 auto numDiff = [&](int i) {
333 const double eps = 1e-6;
334 std::vector<double> p{parametersVec};
335 p[i] = parametersVec[i] - eps;
336 double funcValDown = func(p);
337 p[i] = parametersVec[i] + eps;
338 double funcValUp = func(p);
339 return (funcValUp - funcValDown) / (2 * eps);
340 };
341
342 for (std::size_t i = 0; i < parametersVec.size(); ++i) {
343 std::cout << i << ":" << std::endl;
344 std::cout << " numr : " << numDiff(i) << std::endl;
345 std::cout << " clad : " << gradientVec[i] << std::endl;
346 }
347}
348)";
349}
350
351} // namespace Experimental
352
353} // namespace RooFit
#define g(i)
Definition RSha256.hxx:105
#define oocoutE(o, a)
#define coutE(a)
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
char name[80]
Definition TGX11.cxx:110
#define gInterpreter
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:79
RooFit::OwningPtr< RooArgSet > getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
A class to maintain the context for squashing of RooFit models into code.
std::string assembleCode(std::string const &returnExpr)
Assemble and return the final code with the return expression and global statements.
void addResult(RooAbsArg const *key, std::string const &value)
A function to save an expression that includes/depends on the result of the input node.
void addVecObs(const char *key, int idx)
Since the squashed code represents all observables as a single flattened array, it is important to ke...
std::string const & getResult(RooAbsArg const &arg)
Gets the result for the given node using the node name.
A wrapper class to store a C++ function of type 'double (*)(double*, double*)'.
double(*)(double *, double const *, double const *) Func
std::unique_ptr< RooAbsReal > _absReal
std::string buildCode(RooAbsReal const &head)
void loadParamsAndData(RooAbsArg const *head, RooArgSet const &paramSet, const RooAbsData *data, RooSimultaneous const *simPdf)
std::map< RooFit::Detail::DataKey, ObsInfo > _obsInfos
void(*)(double *, double const *, double const *, double *) Grad
std::string declareFunction(std::string const &funcBody)
RooFuncWrapper(const char *name, const char *title, RooAbsReal &obj, const RooAbsData *data=nullptr, RooSimultaneous const *simPdf=nullptr, bool useEvaluator=false)
Facilitates simultaneous fitting of multiple PDFs to subsets of a given dataset.
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
@ InputArguments