ROOT   Reference Guide
RooBernstein.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * Kyle Cranmer *
7 * *
8 *****************************************************************************/
9
10/** \class RooBernstein
11 \ingroup Roofit
12
13Bernstein basis polynomials are positive-definite in the range [0,1].
14In this implementation, we extend [0,1] to be the range of the parameter.
15There are n+1 Bernstein basis polynomials of degree n:
16\f[
17 B_{i,n}(x) = \begin{pmatrix}n \\\ i \end{pmatrix} x^i \cdot (1-x)^{n-i}
18\f]
19Thus, by providing n coefficients that are positive-definite, there
20is a natural way to have well-behaved polynomial PDFs. For any n, the n+1 polynomials
21'form a partition of unity', i.e., they sum to one for all values of x.
22They can be used as a basis to span the space of polynomials with degree n or less:
23\f[
24 PDF(x, c_0, ..., c_n) = \mathcal{N} \cdot \sum_{i=0}^{n} c_i \cdot B_{i,n}(x).
25\f]
26By giving n+1 coefficients in the constructor, this class constructs the n+1
27polynomials of degree n, and sums them to form an element of the space of polynomials
28of degree n. \f$\mathcal{N} \f$ is a normalisation constant that takes care of the
29cases where the \f$c_i \f$ are not all equal to one.
30
32http://www.idav.ucdavis.edu/education/CAGDNotes/Bernstein-Polynomials.pdf
33**/
34
35#include "RooBernstein.h"
36#include "RooFit.h"
37#include "RooAbsReal.h"
38#include "RooRealVar.h"
39#include "RooArgList.h"
40
41#include "TMath.h"
42
43#include <cmath>
44using namespace std;
45
47
48////////////////////////////////////////////////////////////////////////////////
49
51{
52}
53
54////////////////////////////////////////////////////////////////////////////////
55/// Constructor
56
57RooBernstein::RooBernstein(const char* name, const char* title,
58 RooAbsRealLValue& x, const RooArgList& coefList):
59 RooAbsPdf(name, title),
60 _x("x", "Dependent", this, x),
61 _coefList("coefficients","List of coefficients",this)
62{
63 TIterator* coefIter = coefList.createIterator() ;
64 RooAbsArg* coef ;
65 while((coef = (RooAbsArg*)coefIter->Next())) {
66 if (!dynamic_cast<RooAbsReal*>(coef)) {
67 cout << "RooBernstein::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName()
68 << " is not of type RooAbsReal" << endl ;
69 R__ASSERT(0) ;
70 }
72 }
73 delete coefIter ;
74}
75
76////////////////////////////////////////////////////////////////////////////////
77
78RooBernstein::RooBernstein(const RooBernstein& other, const char* name) :
79 RooAbsPdf(other, name),
80 _x("x", this, other._x),
81 _coefList("coefList",this,other._coefList)
82{
83}
84
85////////////////////////////////////////////////////////////////////////////////
86
87/// Force use of a given normalisation range.
88/// Needed for functions or PDFs (e.g. RooAddPdf) whose shape depends on the choice of normalisation.
89void RooBernstein::selectNormalizationRange(const char* rangeName, Bool_t force)
90{
91 if (rangeName && (force || !_refRangeName.empty())) {
92 _refRangeName = rangeName;
93 }
94}
95
96////////////////////////////////////////////////////////////////////////////////
97
99{
100 double xmax,xmin;
101 std::tie(xmin, xmax) = _x->getRange(_refRangeName.empty() ? nullptr : _refRangeName.c_str());
102 Double_t x = (_x - xmin) / (xmax - xmin); // rescale to [0,1]
103 Int_t degree = _coefList.getSize() - 1; // n+1 polys of degree n
105
106 if(degree == 0) {
107
108 return ((RooAbsReal *)iter.next())->getVal();
109
110 } else if(degree == 1) {
111
112 Double_t a0 = ((RooAbsReal *)iter.next())->getVal(); // c0
113 Double_t a1 = ((RooAbsReal *)iter.next())->getVal() - a0; // c1 - c0
114 return a1 * x + a0;
115
116 } else if(degree == 2) {
117
118 Double_t a0 = ((RooAbsReal *)iter.next())->getVal(); // c0
119 Double_t a1 = 2 * (((RooAbsReal *)iter.next())->getVal() - a0); // 2 * (c1 - c0)
120 Double_t a2 = ((RooAbsReal *)iter.next())->getVal() - a1 - a0; // c0 - 2 * c1 + c2
121 return (a2 * x + a1) * x + a0;
122
123 } else if(degree > 2) {
124
125 Double_t t = x;
126 Double_t s = 1 - x;
127
128 Double_t result = ((RooAbsReal *)iter.next())->getVal() * s;
129 for(Int_t i = 1; i < degree; i++) {
130 result = (result + t * TMath::Binomial(degree, i) * ((RooAbsReal *)iter.next())->getVal()) * s;
131 t *= x;
132 }
133 result += t * ((RooAbsReal *)iter.next())->getVal();
134
135 return result;
136 }
137
138 // in case list of arguments passed is empty
139 return TMath::SignalingNaN();
140}
141
142////////////////////////////////////////////////////////////////////////////////
143
144namespace {
145//Author: Emmanouil Michalainas, CERN 16 AUGUST 2019
146
147void compute( size_t batchSize, double xmax, double xmin,
148 double * __restrict output,
149 const double * __restrict const xData,
150 const RooListProxy& coefList)
151{
152 constexpr size_t block = 128;
153 const int nCoef = coefList.size();
154 const int degree = nCoef-1;
155 double X[block], _1_X[block], powX[block], pow_1_X[block];
156 double *Binomial = new double[nCoef+5];
157 //Binomial stores values c(degree,i) for i in [0..degree]
158
159 Binomial[0] = 1.0;
160 for (int i=1; i<=degree; i++) {
161 Binomial[i] = Binomial[i-1]*(degree-i+1)/i;
162 }
163
164 for (size_t i=0; i<batchSize; i+=block) {
165 const size_t stop = (i+block > batchSize) ? batchSize-i : block;
166
167 //initialization
168 for (size_t j=0; j<stop; j++) {
169 powX[j] = pow_1_X[j] = 1.0;
170 X[j] = (xData[i+j]-xmin) / (xmax-xmin);
171 _1_X[j] = 1-X[j];
172 output[i+j] = 0.0;
173 }
174
175 //raising 1-x to the power of degree
176 for (int k=2; k<=degree; k+=2)
177 for (size_t j=0; j<stop; j++)
178 pow_1_X[j] *= _1_X[j]*_1_X[j];
179
180 if (degree%2 == 1)
181 for (size_t j=0; j<stop; j++)
182 pow_1_X[j] *= _1_X[j];
183
184 //inverting 1-x ---> 1/(1-x)
185 for (size_t j=0; j<stop; j++)
186 _1_X[j] = 1/_1_X[j];
187
188 for (int k=0; k<nCoef; k++) {
189 double coef = static_cast<RooAbsReal&>(coefList[k]).getVal();
190 for (size_t j=0; j<stop; j++) {
191 output[i+j] += coef*Binomial[k]*powX[j]*pow_1_X[j];
192
193 //calculating next power for x and 1-x
194 powX[j] *= X[j];
195 pow_1_X[j] *= _1_X[j];
196 }
197 }
198 }
199 delete[] Binomial;
200}
201};
202
203////////////////////////////////////////////////////////////////////////////////
204
205RooSpan<double> RooBernstein::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
206 auto xData = _x.getValBatch(begin, batchSize);
207 if (xData.empty()) {
208 return {};
209 }
210
211 batchSize = xData.size();
212 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
213 const double xmax = _x.max();
214 const double xmin = _x.min();
215 compute(batchSize, xmax, xmin, output.data(), xData.data(), _coefList);
216 return output;
217}
218
219////////////////////////////////////////////////////////////////////////////////
220
221Int_t RooBernstein::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
222{
223
224 if (matchArgs(allVars, analVars, _x)) return 1;
225 return 0;
226}
227
228////////////////////////////////////////////////////////////////////////////////
229
230Double_t RooBernstein::analyticalIntegral(Int_t code, const char* rangeName) const
231{
232 R__ASSERT(code==1) ;
233
234 double xmax,xmin;
235 std::tie(xmin, xmax) = _x->getRange(_refRangeName.empty() ? nullptr : _refRangeName.c_str());
236 const Double_t xlo = (_x.min(rangeName) - xmin) / (xmax - xmin);
237 const Double_t xhi = (_x.max(rangeName) - xmin) / (xmax - xmin);
238
239 Int_t degree= _coefList.getSize()-1; // n+1 polys of degree n
240 Double_t norm(0) ;
241
243 Double_t temp=0;
244 for (int i=0; i<=degree; ++i){
245 // for each of the i Bernstein basis polynomials
246 // represent it in the 'power basis' (the naive polynomial basis)
247 // where the integral is straight forward.
248 temp = 0;
249 for (int j=i; j<=degree; ++j){ // power basis≈ß
250 temp += pow(-1.,j-i) * TMath::Binomial(degree, j) * TMath::Binomial(j,i) * (TMath::Power(xhi,j+1) - TMath::Power(xlo,j+1)) / (j+1);
251 }
252 temp *= ((RooAbsReal*)iter.next())->getVal(); // include coeff
253 norm += temp; // add this basis's contribution to total
254 }
255
256 norm *= xmax-xmin;
257 return norm;
258}
#define ClassImp(name)
Definition: Rtypes.h:361
#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, const RooArgSet *const normSet=nullptr, Tag_t ownerTag=kUnspecified)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:118
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:73
Int_t getSize() const
RooFIter fwdIterator() const
One-time forward iterator.
Storage_t::size_type size() const
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
std::pair< double, double > getRange(const char *name=0) const
Get low and high bound of the variable.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
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:90
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:450
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
Bernstein basis polynomials are positive-definite in the range [0,1].
Definition: RooBernstein.h:25
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
std::string _refRangeName
Definition: RooBernstein.h:44
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
RooTemplateProxy< RooAbsRealLValue > _x
Definition: RooBernstein.h:42
void selectNormalizationRange(const char *rangeName=0, Bool_t force=kFALSE)
Force use of a given normalisation range.
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooListProxy _coefList
Definition: RooBernstein.h:43
A one-time forward iterator working on RooLinkedList or RooAbsCollection.
RooAbsArg * next()
Return next element or nullptr if at end.
RooListProxy is the concrete proxy for RooArgList objects.
Definition: RooListProxy.h:25
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
double min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
double max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
RooSpan< const typename T::value_type > getValBatch(std::size_t begin, std::size_t batchSize) const
Retrieve a batch of real or category data.
Iterator abstract base class.
Definition: TIterator.h:30
virtual TObject * Next()=0
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Double_t x[n]
Definition: legend1.C:17
constexpr size_t block
Definition: BatchHelpers.h:29
static constexpr double s
static constexpr double degree
Double_t Binomial(Int_t n, Int_t k)
Calculate the binomial coefficient n over k.
Definition: TMath.cxx:2083
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN)
Definition: TMath.h:898
static void output(int code)
Definition: gifencode.c:226