Logo ROOT   master
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 
13 Bernstein basis polynomials are positive-definite in the range [0,1].
14 In this implementation, we extend [0,1] to be the range of the parameter.
15 There 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]
19 Thus, by providing n coefficients that are positive-definite, there
20 is 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.
22 They 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]
26 By giving n+1 coefficients in the constructor, this class constructs the n+1
27 polynomials of degree n, and sums them to form an element of the space of polynomials
28 of degree n. \f$ \mathcal{N} \f$ is a normalisation constant that takes care of the
29 cases where the \f$ c_i \f$ are not all equal to one.
30 
31 See also
32 http://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>
44 using namespace std;
45 
47 
48 ////////////////////////////////////////////////////////////////////////////////
49 
51 {
52 }
53 
54 ////////////////////////////////////////////////////////////////////////////////
55 /// Constructor
56 
57 RooBernstein::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 
78 RooBernstein::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 
132 namespace {
133 //Author: Emmanouil Michalainas, CERN 16 AUGUST 2019
134 
135 void 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 
193 RooSpan<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 
210 Int_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 
222 Double_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 
229  RooFIter iter = _coefList.fwdIterator() ;
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 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
RooRealProxy _x
Definition: RooBernstein.h:39
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
float xmin
Definition: THbookFile.cxx:93
Storage_t::size_type size() const
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
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
RooListProxy _coefList
Definition: RooBernstein.h:40
#define R__ASSERT(e)
Definition: TError.h:96
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
No analytical calculation available (yet) of integrals over subranges.
STL namespace.
Iterator abstract base class.
Definition: TIterator.h:30
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral. ...
Double_t x[n]
Definition: legend1.C:17
double pow(double, double)
RooSpan< const typename T::value_type > getValBatch(std::size_t begin, std::size_t batchSize) const
Retrieve a batch of real or category data.
Definition: RooProxy.h:240
TIterator * createIterator(Bool_t dir=kIterForward) const R__SUGGEST_ALTERNATIVE("begin()
TIterator-style iteration over contained elements.
static constexpr double s
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
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
double min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
Definition: RooProxy.h:262
static constexpr double degree
Int_t getSize() const
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:89
Double_t Binomial(Int_t n, Int_t k)
float xmax
Definition: THbookFile.cxx:93
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN) ...
Definition: TMath.h:898
RooListProxy is the concrete proxy for RooArgList objects.
Definition: RooListProxy.h:25
RooAbsArg * next()
Return next element or nullptr if at end.
#define ClassImp(name)
Definition: Rtypes.h:365
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
virtual TObject * Next()=0
constexpr size_t block
Definition: BatchHelpers.h:29
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:70
static void output(int code)
Definition: gifencode.c:226
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:449
RooFIter fwdIterator() const R__SUGGEST_ALTERNATIVE("begin()
One-time forward iterator.
char name[80]
Definition: TGX11.cxx:109
A one-time forward iterator working on RooLinkedList or RooAbsCollection.
double max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
Definition: RooProxy.h:264