Logo 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
31See also
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 RooAbsReal& 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 }
71 _coefList.add(*coef) ;
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
88{
89 Double_t xmin = _x.min();
90 Double_t x = (_x - xmin) / (_x.max() - xmin); // rescale to [0,1]
91 Int_t degree = _coefList.getSize() - 1; // n+1 polys of degree n
93
94 if(degree == 0) {
95
96 return ((RooAbsReal *)iter.next())->getVal();
97
98 } else if(degree == 1) {
99
100 Double_t a0 = ((RooAbsReal *)iter.next())->getVal(); // c0
101 Double_t a1 = ((RooAbsReal *)iter.next())->getVal() - a0; // c1 - c0
102 return a1 * x + a0;
103
104 } else if(degree == 2) {
105
106 Double_t a0 = ((RooAbsReal *)iter.next())->getVal(); // c0
107 Double_t a1 = 2 * (((RooAbsReal *)iter.next())->getVal() - a0); // 2 * (c1 - c0)
108 Double_t a2 = ((RooAbsReal *)iter.next())->getVal() - a1 - a0; // c0 - 2 * c1 + c2
109 return (a2 * x + a1) * x + a0;
110
111 } else if(degree > 2) {
112
113 Double_t t = x;
114 Double_t s = 1 - x;
115
116 Double_t result = ((RooAbsReal *)iter.next())->getVal() * s;
117 for(Int_t i = 1; i < degree; i++) {
118 result = (result + t * TMath::Binomial(degree, i) * ((RooAbsReal *)iter.next())->getVal()) * s;
119 t *= x;
120 }
121 result += t * ((RooAbsReal *)iter.next())->getVal();
122
123 return result;
124 }
125
126 // in case list of arguments passed is empty
127 return TMath::SignalingNaN();
128}
129
130////////////////////////////////////////////////////////////////////////////////
131
132namespace {
133//Author: Emmanouil Michalainas, CERN 16 AUGUST 2019
134
135void compute( size_t batchSize, double xmax, double xmin,
136 double * __restrict output,
137 const double * __restrict const xData,
138 const RooListProxy& coefList)
139{
140 constexpr size_t block = 128;
141 const int nCoef = coefList.size();
142 const int degree = nCoef-1;
143 double X[block], _1_X[block], powX[block], pow_1_X[block];
144 double *Binomial = new double[nCoef+5];
145 //Binomial stores values c(degree,i) for i in [0..degree]
146
147 Binomial[0] = 1.0;
148 for (int i=1; i<=degree; i++) {
149 Binomial[i] = Binomial[i-1]*(degree-i+1)/i;
150 }
151
152 for (size_t i=0; i<batchSize; i+=block) {
153 const size_t stop = (i+block > batchSize) ? batchSize-i : block;
154
155 //initialization
156 for (size_t j=0; j<stop; j++) {
157 powX[j] = pow_1_X[j] = 1.0;
158 X[j] = (xData[i+j]-xmin) / (xmax-xmin);
159 _1_X[j] = 1-X[j];
160 output[i+j] = 0.0;
161 }
162
163 //raising 1-x to the power of degree
164 for (int k=2; k<=degree; k+=2)
165 for (size_t j=0; j<stop; j++)
166 pow_1_X[j] *= _1_X[j]*_1_X[j];
167
168 if (degree%2 == 1)
169 for (size_t j=0; j<stop; j++)
170 pow_1_X[j] *= _1_X[j];
171
172 //inverting 1-x ---> 1/(1-x)
173 for (size_t j=0; j<stop; j++)
174 _1_X[j] = 1/_1_X[j];
175
176 for (int k=0; k<nCoef; k++) {
177 double coef = static_cast<RooAbsReal&>(coefList[k]).getVal();
178 for (size_t j=0; j<stop; j++) {
179 output[i+j] += coef*Binomial[k]*powX[j]*pow_1_X[j];
180
181 //calculating next power for x and 1-x
182 powX[j] *= X[j];
183 pow_1_X[j] *= _1_X[j];
184 }
185 }
186 }
187 delete[] Binomial;
188}
189};
190
191////////////////////////////////////////////////////////////////////////////////
192
193RooSpan<double> RooBernstein::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
194 auto xData = _x.getValBatch(begin, batchSize);
195 if (xData.empty()) {
196 return {};
197 }
198
199 batchSize = xData.size();
200 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
201 const double xmax = _x.max();
202 const double xmin = _x.min();
203 compute(batchSize, xmax, xmin, output.data(), xData.data(), _coefList);
204 return output;
205}
206
207////////////////////////////////////////////////////////////////////////////////
208/// No analytical calculation available (yet) of integrals over subranges
209
210Int_t RooBernstein::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
211{
212 if (rangeName && strlen(rangeName)) {
213 return 0 ;
214 }
215
216 if (matchArgs(allVars, analVars, _x)) return 1;
217 return 0;
218}
219
220////////////////////////////////////////////////////////////////////////////////
221
222Double_t RooBernstein::analyticalIntegral(Int_t code, const char* rangeName) const
223{
224 R__ASSERT(code==1) ;
225 Double_t xmin = _x.min(rangeName); Double_t xmax = _x.max(rangeName);
226 Int_t degree= _coefList.getSize()-1; // n+1 polys of degree n
227 Double_t norm(0) ;
228
230 Double_t temp=0;
231 for (int i=0; i<=degree; ++i){
232 // for each of the i Bernstein basis polynomials
233 // represent it in the 'power basis' (the naive polynomial basis)
234 // where the integral is straight forward.
235 temp = 0;
236 for (int j=i; j<=degree; ++j){ // power basisŧ
237 temp += pow(-1.,j-i) * TMath::Binomial(degree, j) * TMath::Binomial(j,i) / (j+1);
238 }
239 temp *= ((RooAbsReal*)iter.next())->getVal(); // include coeff
240 norm += temp; // add this basis's contribution to total
241 }
242
243 norm *= xmax-xmin;
244 return norm;
245}
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
Storage_t::size_type size() const
TIterator * createIterator(Bool_t dir=kIterForward) const R__SUGGEST_ALTERNATIVE("begin()
TIterator-style iteration over contained elements.
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
Bernstein basis polynomials are positive-definite in the range [0,1].
Definition: RooBernstein.h:23
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
No analytical calculation available (yet) of integrals over subranges.
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
RooRealProxy _x
Definition: RooBernstein.h:39
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooListProxy _coefList
Definition: RooBernstein.h:40
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)
Reimplementation of standard RooArgList::add()
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
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.
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:2090
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