Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooPolyFunc.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2021: *
10 * CERN, Switzerland *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/** \class RooPolyFunc
18 \ingroup Roofit
19
20RooPolyFunc implements a polynomial function in multi-variables.
21The polynomial coefficients are implemented as doubles and are not
22part of the RooFit computation graph.
23**/
24
25#include "RooPolyFunc.h"
26
27#include "RooAbsReal.h"
28#include "RooArgList.h"
29#include "RooArgSet.h"
30#include "RooDerivative.h"
31#include "RooMsgService.h"
32#include "RooRealVar.h"
33
34#include <utility>
35
36using namespace std;
37using namespace RooFit;
38
40
41////////////////////////////////////////////////////////////////////////////////
42/// coverity[UNINIT_CTOR]
43
44void RooPolyFunc::addTerm(double coefficient)
45{
46 int n_terms = _terms.size();
47 std::string coeff_name = Form("%s_c%d", GetName(), n_terms);
48 std::string term_name = Form("%s_t%d", GetName(), n_terms);
49 auto termList = std::make_unique<RooListProxy>(term_name.c_str(), term_name.c_str(), this);
50 auto coeff = std::make_unique<RooRealVar>(coeff_name.c_str(), coeff_name.c_str(), coefficient);
51 RooArgList exponents{};
52
53 for (const auto &var : _vars) {
54 std::string exponent_name = Form("%s_%s^%d", GetName(), var->GetName(), 0);
55 auto exponent = std::make_unique<RooRealVar>(exponent_name.c_str(), exponent_name.c_str(), 0);
56 exponents.addOwned(std::move(exponent));
57 }
58
59 termList->addOwned(std::move(exponents));
60 termList->addOwned(std::move(coeff));
61 _terms.push_back(std::move(termList));
62}
63
64void RooPolyFunc::addTerm(double coefficient, const RooAbsReal &var1, int exp1)
65{
66 int n_terms = _terms.size();
67 std::string coeff_name = Form("%s_c%d", GetName(), n_terms);
68 std::string term_name = Form("%s_t%d", GetName(), n_terms);
69 auto termList = std::make_unique<RooListProxy>(term_name.c_str(), term_name.c_str(), this);
70 auto coeff = std::make_unique<RooRealVar>(coeff_name.c_str(), coeff_name.c_str(), coefficient);
71 RooArgList exponents{};
72
73 // linear iterate over all the variables, create var1^exp1 ..vark^0
74 for (const auto &var : _vars) {
75 int exp = 0;
76 if (strcmp(var1.GetName(), var->GetName()) == 0)
77 exp += exp1;
78 std::string exponent_name = Form("%s_%s^%d", GetName(), var->GetName(), exp);
79 auto exponent = std::make_unique<RooRealVar>(exponent_name.c_str(), exponent_name.c_str(), exp);
80 exponents.addOwned(std::move(exponent));
81 }
82
83 termList->addOwned(std::move(exponents));
84 termList->addOwned(std::move(coeff));
85 _terms.push_back(std::move(termList));
86}
87
88void RooPolyFunc::addTerm(double coefficient, const RooAbsReal &var1, int exp1, const RooAbsReal &var2, int exp2)
89{
90
91 int n_terms = _terms.size();
92 std::string coeff_name = Form("%s_c%d", GetName(), n_terms);
93 std::string term_name = Form("%s_t%d", GetName(), n_terms);
94 auto termList = std::make_unique<RooListProxy>(term_name.c_str(), term_name.c_str(), this);
95 auto coeff = std::make_unique<RooRealVar>(coeff_name.c_str(), coeff_name.c_str(), coefficient);
96 RooArgList exponents{};
97
98 for (const auto &var : _vars) {
99 int exp = 0;
100 if (strcmp(var1.GetName(), var->GetName()) == 0)
101 exp += exp1;
102 if (strcmp(var2.GetName(), var->GetName()) == 0)
103 exp += exp2;
104 std::string exponent_name = Form("%s_%s^%d", GetName(), var->GetName(), exp);
105 auto exponent = std::make_unique<RooRealVar>(exponent_name.c_str(), exponent_name.c_str(), exp);
106 exponents.addOwned(std::move(exponent));
107 }
108 termList->addOwned(std::move(exponents));
109 termList->addOwned(std::move(coeff));
110 _terms.push_back(std::move(termList));
111}
112
113void RooPolyFunc::addTerm(double coefficient, const RooAbsCollection &exponents)
114{
115 if (exponents.size() != _vars.size()) {
116 coutE(InputArguments) << "RooPolyFunc::addTerm(" << GetName() << ") WARNING: number of exponents ("
117 << exponents.size() << ") provided do not match the number of variables (" << _vars.size()
118 << ")" << endl;
119 }
120 int n_terms = _terms.size();
121 std::string coeff_name = Form("%s_c%d", GetName(), n_terms);
122 std::string term_name = Form("%s_t%d", GetName(), n_terms);
123 auto termList = std::make_unique<RooListProxy>(term_name.c_str(), term_name.c_str(), this);
124 auto coeff = std::make_unique<RooRealVar>(coeff_name.c_str(), coeff_name.c_str(), coefficient);
125 termList->addOwned(exponents);
126 termList->addOwned(std::move(coeff));
127 _terms.push_back(std::move(termList));
128}
129
130////////////////////////////////////////////////////////////////////////////////
131/// Default constructor
132
134
135////////////////////////////////////////////////////////////////////////////////
136/// Parameterised constructor
137
138RooPolyFunc::RooPolyFunc(const char *name, const char *title, const RooAbsCollection &vars)
139 : RooAbsReal(name, title), _vars("x", "list of dependent variables", this)
140{
141 for (const auto &var : vars) {
142 if (!dynamic_cast<RooAbsReal *>(var)) {
143 std::stringstream ss;
144 ss << "RooPolyFunc::ctor(" << GetName() << ") ERROR: coefficient " << var->GetName()
145 << " is not of type RooAbsReal";
146 const std::string errorMsg = ss.str();
147 coutE(InputArguments) << errorMsg << std::endl;
148 throw std::runtime_error(errorMsg);
149 }
150 _vars.add(*var);
151 }
152}
153
154////////////////////////////////////////////////////////////////////////////////
155/// Copy constructor
156
157RooPolyFunc::RooPolyFunc(const RooPolyFunc &other, const char *name)
158 : RooAbsReal(other, name), _vars("vars", this, other._vars)
159{
160 for (auto const &term : other._terms) {
161 _terms.emplace_back(std::make_unique<RooListProxy>(term->GetName(), this, *term));
162 }
163}
164
165////////////////////////////////////////////////////////////////////////////////
166/// Return to RooPolyFunc as a string
167
168std::string RooPolyFunc::asString() const
169{
170 std::stringstream ss;
171 bool first = true;
172 for (const auto &term : _terms) {
173 size_t n_vars = term->size() - 1;
174 auto coef = dynamic_cast<RooRealVar *>(term->at(n_vars));
175 if (coef->getVal() > 0 && !first)
176 ss << "+";
177 ss << coef->getVal();
178 first = true;
179 for (size_t i_var = 0; i_var < n_vars; ++i_var) {
180 auto var = dynamic_cast<RooRealVar *>(_vars.at(i_var));
181 auto exp = dynamic_cast<RooRealVar *>(term->at(i_var));
182 if (exp->getVal() == 0)
183 continue;
184 if (first) {
185 ss << " * (";
186 } else {
187 ss << "*";
188 }
189 ss << "pow(" << var->GetName() << "," << exp->getVal() << ")";
190 first = false;
191 }
192 if (!first)
193 ss << ")";
194 }
195 return ss.str();
196}
197
198////////////////////////////////////////////////////////////////////////////////
199/// Evaluate value of Polynomial.
201{
202 // Calculate and return value of polynomial
203 double poly_sum(0.0);
204 for (const auto &term : _terms) {
205 double poly_term(1.0);
206 size_t n_vars = term->size() - 1;
207 for (size_t i_var = 0; i_var < n_vars; ++i_var) {
208 auto var = dynamic_cast<RooRealVar *>(_vars.at(i_var));
209 auto exp = dynamic_cast<RooRealVar *>(term->at(i_var));
210 poly_term *= pow(var->getVal(), exp->getVal());
211 }
212 auto coef = dynamic_cast<RooRealVar *>(term->at(n_vars));
213 poly_sum += coef->getVal() * poly_term;
214 }
215 return poly_sum;
216}
217
218void setCoordinates(const RooAbsCollection &observables, std::vector<double> const &observableValues)
219// set all observables to expansion co-ordinate
220{
221 std::size_t iObs = 0;
222 for (auto *var : static_range_cast<RooRealVar *>(observables)) {
223 var->setVal(observableValues[iObs++]);
224 }
225}
226
227void fixObservables(const RooAbsCollection &observables)
228{
229 for (auto *var : static_range_cast<RooRealVar *>(observables)) {
230 var->setConstant(true);
231 }
232}
233
234//////////////////////////////////////////////////////////////////////////////
235/// Taylor expanding given function in terms of observables around
236/// observableValues. Supports expansions upto order 2.
237/// \param[in] func Function of variables that is taylor expanded.
238/// \param[in] observables Set of variables to perform the expansion.
239/// It's type is RooArgList to ensure that it is always ordered the
240/// same as the observableValues vector. However, duplicate
241/// observables are still not allowed.
242/// \param[in] order Order of the expansion (0,1,2 supported).
243/// \param[in] observableValues Coordinates around which expansion is
244/// performed. If empty, the nominal observable values are taken, if
245/// the size matches the size of the observables RooArgSet, the
246/// values are mapped to the observables in matching order. If it
247/// contains only one element, the same single value is used for all
248/// observables.
249/// \param[in] eps1 Precision for first derivative and second derivative.
250/// \param[in] eps2 Precision for second partial derivative of cross-derivative.
251std::unique_ptr<RooPolyFunc>
252RooPolyFunc::taylorExpand(const char *name, const char *title, RooAbsReal &func, const RooArgList &observables,
253 int order, std::vector<double> const &observableValues, double eps1, double eps2)
254{
255 // Create the taylor expansion polynomial
256 auto taylorPoly = std::make_unique<RooPolyFunc>(name, title, observables);
257
258 // Verify that there are no duplicate observables
259 {
260 RooArgSet obsSet;
261 for (RooAbsArg *obs : observables) {
262 obsSet.add(*obs, /*silent*/ true); // we can be silent now, the error will come later
263 }
264 if (obsSet.size() != observables.size()) {
265 std::stringstream errorMsgStream;
266 errorMsgStream << "RooPolyFunc::taylorExpand(" << name << ") ERROR: duplicate input observables!";
267 const auto errorMsg = errorMsgStream.str();
268 oocoutE(taylorPoly.get(), InputArguments) << errorMsg << std::endl;
269 throw std::invalid_argument(errorMsg);
270 }
271 }
272
273 // Figure out the observable values around which to exapnd
274 std::vector<double> obsValues;
275 if (observableValues.empty()) {
276 obsValues.reserve(observables.size());
277 for (auto *var : static_range_cast<RooRealVar *>(observables)) {
278 obsValues.push_back(var->getVal());
279 }
280 } else if (observableValues.size() == 1) {
281 obsValues.resize(observables.size());
282 std::fill(obsValues.begin(), obsValues.end(), observableValues[0]);
283 } else if (observableValues.size() == observables.size()) {
284 obsValues = observableValues;
285 } else {
286 std::stringstream errorMsgStream;
287 errorMsgStream << "RooPolyFunc::taylorExpand(" << name
288 << ") ERROR: observableValues must be empty, contain one element, or match observables.size()!";
289 const auto errorMsg = errorMsgStream.str();
290 oocoutE(taylorPoly.get(), InputArguments) << errorMsg << std::endl;
291 throw std::invalid_argument(errorMsg);
292 }
293
294 // taylor expansion can be performed only for order 0, 1, 2 currently
295 if (order >= 3 || order <= 0) {
296 std::stringstream errorMsgStream;
297 errorMsgStream << "RooPolyFunc::taylorExpand(" << name << ") ERROR: order must be 0, 1, or 2";
298 const auto errorMsg = errorMsgStream.str();
299 oocoutE(taylorPoly.get(), InputArguments) << errorMsg << std::endl;
300 throw std::invalid_argument(errorMsg);
301 }
302
303 setCoordinates(observables, obsValues);
304
305 // estimate taylor expansion polynomial for different orders
306 // f(x) = f(x=x0)
307 // + \sum_i + \frac{df}{dx_i}|_{x_i=x_i^0}(x - x_i^0)
308 // + \sum_{i,j} 0.5 * \frac{d^f}{dx_i x_j}
309 // ( x_i - x_i^0 )( x_j - x_j^0 )
310 // note in the polynomial the () brackets are expanded out
311 for (int i_order = 0; i_order <= order; ++i_order) {
312 switch (i_order) {
313 case 0: {
314 taylorPoly->addTerm(func.getVal());
315 break;
316 }
317 case 1: {
318 for (auto *var : static_range_cast<RooRealVar *>(observables)) {
319 double var1_val = var->getVal();
320 auto deriv = func.derivative(*var, 1, eps2);
321 double deriv_val = deriv->getVal();
322 setCoordinates(observables, obsValues);
323 taylorPoly->addTerm(deriv_val, *var, 1);
324 if (var1_val != 0.0) {
325 taylorPoly->addTerm(deriv_val * var1_val * -1.0);
326 }
327 }
328 break;
329 }
330 case 2: {
331 for (auto *var1 : static_range_cast<RooRealVar *>(observables)) {
332 double var1_val = var1->getVal();
333 auto deriv1 = func.derivative(*var1, 1, eps1);
334 for (auto *var2 : static_range_cast<RooRealVar *>(observables)) {
335 double var2_val = var2->getVal();
336 double deriv_val = 0.0;
337 if (strcmp(var1->GetName(), var2->GetName()) == 0) {
338 auto deriv2 = func.derivative(*var2, 2, eps2);
339 deriv_val = 0.5 * deriv2->getVal();
340 setCoordinates(observables, obsValues);
341 } else {
342 auto deriv2 = deriv1->derivative(*var2, 1, eps2);
343 deriv_val = 0.5 * deriv2->getVal();
344 setCoordinates(observables, obsValues);
345 }
346 taylorPoly->addTerm(deriv_val, *var1, 1, *var2, 1);
347 if (var1_val != 0.0 || var2_val != 0.0) {
348 taylorPoly->addTerm(deriv_val * var1_val * var2_val);
349 taylorPoly->addTerm(deriv_val * var2_val * -1.0, *var1, 1);
350 taylorPoly->addTerm(deriv_val * var1_val * -1.0, *var2, 1);
351 }
352 }
353 }
354 break;
355 }
356 }
357 }
358 return taylorPoly;
359}
#define oocoutE(o, a)
#define coutE(a)
void fixObservables(const RooAbsCollection &observables)
void setCoordinates(const RooAbsCollection &observables, std::vector< double > const &observableValues)
#define ClassImp(name)
Definition Rtypes.h:377
char name[80]
Definition TGX11.cxx:110
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2468
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:79
Abstract container object that can hold multiple RooAbsArg objects.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
RooDerivative * derivative(RooRealVar &obs, Int_t order=1, double eps=0.001)
Return function representing first, second or third order derivative of this function.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
RooPolyFunc implements a polynomial function in multi-variables.
Definition RooPolyFunc.h:28
RooListProxy _vars
Definition RooPolyFunc.h:62
double evaluate() const override
Evaluation.
RooPolyFunc()
Default constructor.
void addTerm(double coefficient)
coverity[UNINIT_CTOR]
static std::unique_ptr< RooPolyFunc > taylorExpand(const char *name, const char *title, RooAbsReal &func, const RooArgList &observables, int order=1, std::vector< double > const &observableValues={}, double eps1=1e-6, double eps2=1e-3)
Taylor expanding given function in terms of observables around observableValues.
std::string asString() const
Return to RooPolyFunc as a string.
std::vector< std::unique_ptr< RooListProxy > > _terms
Definition RooPolyFunc.h:63
Variable that can be changed from the outside.
Definition RooRealVar.h:37
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
@ InputArguments
Definition first.py:1