ROOT   Reference Guide
RooPolynomial.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) 2000-2005, Regents of the University of California *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
15 *****************************************************************************/
16
17/** \class RooPolynomial
18 \ingroup Roofit
19
20RooPolynomial implements a polynomial p.d.f of the form
21\f[ f(x) = \mathcal{N} \cdot \sum_{i} a_{i} * x^i \f]
22By default, the coefficient \f$a_0 \f$ is chosen to be 1, as polynomial
23probability density functions have one degree of freedom
24less than polynomial functions due to the normalisation condition. \f$\mathcal{N} \f$
25is a normalisation constant that is automatically calculated when the polynomial is used
26in computations.
27
28The sum can be truncated at the low end. See the main constructor
29RooPolynomial::RooPolynomial(const char*, const char*, RooAbsReal&, const RooArgList&, Int_t)
30**/
31
32#include "RooPolynomial.h"
33#include "RooAbsReal.h"
34#include "RooArgList.h"
35#include "RooMsgService.h"
36#include "BatchHelpers.h"
37
38#include "TError.h"
39
40#include <cmath>
41#include <cassert>
42#include <vector>
43using namespace std;
44
46
47////////////////////////////////////////////////////////////////////////////////
48/// coverity[UNINIT_CTOR]
49
51{
52}
53
54////////////////////////////////////////////////////////////////////////////////
55/// Create a polynomial in the variable x.
56/// \param[in] name Name of the PDF
57/// \param[in] title Title for plotting the PDF
58/// \param[in] x The variable of the polynomial
59/// \param[in] coefList The coefficients \f$a_i \f$
60/// \param[in] lowestOrder [optional] Truncate the sum such that it skips the lower orders:
61/// \f[
62/// 1. + \sum_{i=0}^{\mathrm{coefList.size()}} a_{i} * x^{(i + \mathrm{lowestOrder})}
63/// \f]
64///
65/// This means that
66/// \code{.cpp}
67/// RooPolynomial pol("pol", "pol", x, RooArgList(a, b), lowestOrder = 2)
68/// \endcode
69/// computes
70/// \f[
71/// \mathrm{pol}(x) = 1 * x^0 + (0 * x^{\ldots}) + a * x^2 + b * x^3.
72/// \f]
73
74
75RooPolynomial::RooPolynomial(const char* name, const char* title,
76 RooAbsReal& x, const RooArgList& coefList, Int_t lowestOrder) :
77 RooAbsPdf(name, title),
78 _x("x", "Dependent", this, x),
79 _coefList("coefList","List of coefficients",this),
80 _lowestOrder(lowestOrder)
81{
82 // Check lowest order
83 if (_lowestOrder<0) {
84 coutE(InputArguments) << "RooPolynomial::ctor(" << GetName()
85 << ") WARNING: lowestOrder must be >=0, setting value to 0" << endl ;
86 _lowestOrder=0 ;
87 }
88
89 RooFIter coefIter = coefList.fwdIterator() ;
90 RooAbsArg* coef ;
91 while((coef = (RooAbsArg*)coefIter.next())) {
92 if (!dynamic_cast<RooAbsReal*>(coef)) {
93 coutE(InputArguments) << "RooPolynomial::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName()
94 << " is not of type RooAbsReal" << endl ;
95 R__ASSERT(0) ;
96 }
98 }
99}
100
101////////////////////////////////////////////////////////////////////////////////
102
103RooPolynomial::RooPolynomial(const char* name, const char* title,
104 RooAbsReal& x) :
105 RooAbsPdf(name, title),
106 _x("x", "Dependent", this, x),
107 _coefList("coefList","List of coefficients",this),
108 _lowestOrder(1)
109{ }
110
111////////////////////////////////////////////////////////////////////////////////
112/// Copy constructor
113
115 RooAbsPdf(other, name),
116 _x("x", this, other._x),
117 _coefList("coefList",this,other._coefList),
118 _lowestOrder(other._lowestOrder)
119{ }
120
121////////////////////////////////////////////////////////////////////////////////
122/// Destructor
123
125{ }
126
127////////////////////////////////////////////////////////////////////////////////
128
130{
131 // Calculate and return value of polynomial
132
133 const unsigned sz = _coefList.getSize();
134 const int lowestOrder = _lowestOrder;
135 if (!sz) return lowestOrder ? 1. : 0.;
136 _wksp.clear();
137 _wksp.reserve(sz);
138 {
139 const RooArgSet* nset = _coefList.nset();
141 RooAbsReal* c;
142 while ((c = (RooAbsReal*) it.next())) _wksp.push_back(c->getVal(nset));
143 }
144 const Double_t x = _x;
145 Double_t retVal = _wksp[sz - 1];
146 for (unsigned i = sz - 1; i--; ) retVal = _wksp[i] + x * retVal;
147 return retVal * std::pow(x, lowestOrder) + (lowestOrder ? 1.0 : 0.0);
148}
149
150////////////////////////////////////////////////////////////////////////////////
151
152namespace {
153//Author: Emmanouil Michalainas, CERN 15 AUGUST 2019
154
155void compute( size_t batchSize, const int lowestOrder,
156 double * __restrict output,
157 const double * __restrict const X,
159{
160 const int nCoef = coefList.size();
161 if (nCoef==0 && lowestOrder==0) {
162 for (size_t i=0; i<batchSize; i++) {
163 output[i] = 0.0;
164 }
165 }
166 else if (nCoef==0 && lowestOrder>0) {
167 for (size_t i=0; i<batchSize; i++) {
168 output[i] = 1.0;
169 }
170 } else {
171 for (size_t i=0; i<batchSize; i++) {
172 output[i] = coefList[nCoef-1][i];
173 }
174 }
175 if (nCoef == 0) return;
176
177 /* Indexes are in range 0..nCoef-1 but coefList[nCoef-1]
178 * has already been processed. In order to traverse the list,
179 * with step of 2 we have to start at index nCoef-3 and use
180 * coefList[k+1] and coefList[k]
181 */
182 for (int k=nCoef-3; k>=0; k-=2) {
183 for (size_t i=0; i<batchSize; i++) {
184 double coef1 = coefList[k+1][i];
185 double coef2 = coefList[ k ][i];
186 output[i] = X[i]*(output[i]*X[i] + coef1) + coef2;
187 }
188 }
189 // If nCoef is odd, then the coefList[0] didn't get processed
190 if (nCoef%2 == 0) {
191 for (size_t i=0; i<batchSize; i++) {
192 output[i] = output[i]*X[i] + coefList[0][i];
193 }
194 }
195 //Increase the order of the polynomial, first by myltiplying with X[i]^2
196 if (lowestOrder == 0) return;
197 for (int k=2; k<=lowestOrder; k+=2) {
198 for (size_t i=0; i<batchSize; i++) {
199 output[i] *= X[i]*X[i];
200 }
201 }
202 const bool isOdd = lowestOrder%2;
203 for (size_t i=0; i<batchSize; i++) {
204 if (isOdd) output[i] *= X[i];
205 output[i] += 1.0;
206 }
207}
208};
209
210////////////////////////////////////////////////////////////////////////////////
211
212RooSpan<double> RooPolynomial::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
213 RooSpan<const double> xData = _x.getValBatch(begin, batchSize);
214 batchSize = xData.size();
215 if (xData.empty()) {
216 return {};
217 }
218
219 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
220 const int nCoef = _coefList.getSize();
221 const RooArgSet* normSet = _coefList.nset();
223 for (int i=0; i<nCoef; i++) {
224 auto val = static_cast<RooAbsReal&>(_coefList[i]).getVal(normSet);
225 auto valBatch = static_cast<RooAbsReal&>(_coefList[i]).getValBatch(begin, batchSize, normSet);
226 coefList.emplace_back(val, valBatch);
227 }
228
229 compute(batchSize, _lowestOrder, output.data(), xData.data(), coefList);
230
231 return output;
232}
233
234////////////////////////////////////////////////////////////////////////////////
235/// Advertise to RooFit that this function can be analytically integrated.
236Int_t RooPolynomial::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
237{
238 if (matchArgs(allVars, analVars, _x)) return 1;
239 return 0;
240}
241
242////////////////////////////////////////////////////////////////////////////////
243/// Do the analytical integral according to the code that was returned by getAnalyticalIntegral().
244Double_t RooPolynomial::analyticalIntegral(Int_t code, const char* rangeName) const
245{
246 R__ASSERT(code==1) ;
247
248 const Double_t xmin = _x.min(rangeName), xmax = _x.max(rangeName);
249 const int lowestOrder = _lowestOrder;
250 const unsigned sz = _coefList.getSize();
251 if (!sz) return xmax - xmin;
252 _wksp.clear();
253 _wksp.reserve(sz);
254 {
255 const RooArgSet* nset = _coefList.nset();
257 unsigned i = 1 + lowestOrder;
258 RooAbsReal* c;
259 while ((c = (RooAbsReal*) it.next())) {
260 _wksp.push_back(c->getVal(nset) / Double_t(i));
261 ++i;
262 }
263 }
264 Double_t min = _wksp[sz - 1], max = _wksp[sz - 1];
265 for (unsigned i = sz - 1; i--; )
266 min = _wksp[i] + xmin * min, max = _wksp[i] + xmax * max;
267 return max * std::pow(xmax, 1 + lowestOrder) - min * std::pow(xmin, 1 + lowestOrder) +
268 (lowestOrder ? (xmax - xmin) : 0.);
269}
#define c(i)
Definition: RSha256.hxx:101
#define coutE(a)
Definition: RooMsgService.h:34
int Int_t
Definition: RtypesCore.h:41
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
#define R__ASSERT(e)
Definition: TError.h:96
char name[80]
Definition: TGX11.cxx:109
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
double pow(double, double)
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:58
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:71
RooFIter fwdIterator() const R__SUGGEST_ALTERNATIVE("begin()
One-time forward iterator.
Int_t getSize() const
virtual RooSpan< const double > getValBatch(std::size_t begin, std::size_t batchSize, const RooArgSet *normSet=nullptr) const
Compute batch of values for given range, and normalise by integrating over the observables in nset.
Definition: RooAbsPdf.cxx:348
const RooArgSet * nset() const
Definition: RooAbsProxy.h:46
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:87
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:447
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
A one-time forward iterator working on RooLinkedList or RooAbsCollection.
RooAbsArg * next()
Return next element or nullptr if at end.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
RooPolynomial implements a polynomial p.d.f of the form.
Definition: RooPolynomial.h:28
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Advertise to RooFit that this function can be analytically integrated.
RooRealProxy _x
Definition: RooPolynomial.h:45
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Do the analytical integral according to the code that was returned by getAnalyticalIntegral().
Int_t _lowestOrder
Definition: RooPolynomial.h:47
RooPolynomial()
coverity[UNINIT_CTOR]
Double_t evaluate() const
do not persist
RooListProxy _coefList
Definition: RooPolynomial.h:46
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
std::vector< Double_t > _wksp
Definition: RooPolynomial.h:49
virtual ~RooPolynomial()
Destructor.
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
constexpr std::span< T >::pointer data() const
Definition: RooSpan.h:115
constexpr std::span< T >::index_type size() const noexcept
Definition: RooSpan.h:123
constexpr bool empty() const noexcept
Definition: RooSpan.h:127
Double_t min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
RooSpan< const double > getValBatch(std::size_t begin, std::size_t batchSize) const
Double_t max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Double_t x[n]
Definition: legend1.C:17
@ InputArguments
Definition: RooGlobalFunc.h:68
static void output(int code)
Definition: gifencode.c:226