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#include "RooBatchCompute.h"
41
42#include "TMath.h"
43
44#include <cmath>
45using namespace std;
46
48
49////////////////////////////////////////////////////////////////////////////////
50
52{
53}
54
55////////////////////////////////////////////////////////////////////////////////
56/// Constructor
57
58RooBernstein::RooBernstein(const char* name, const char* title,
59 RooAbsRealLValue& x, const RooArgList& coefList):
60 RooAbsPdf(name, title),
61 _x("x", "Dependent", this, x),
62 _coefList("coefficients","List of coefficients",this)
63{
64 TIterator* coefIter = coefList.createIterator() ;
65 RooAbsArg* coef ;
66 while((coef = (RooAbsArg*)coefIter->Next())) {
67 if (!dynamic_cast<RooAbsReal*>(coef)) {
68 cout << "RooBernstein::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName()
69 << " is not of type RooAbsReal" << endl ;
70 R__ASSERT(0) ;
71 }
73 }
74 delete coefIter ;
75}
76
77////////////////////////////////////////////////////////////////////////////////
78
79RooBernstein::RooBernstein(const RooBernstein& other, const char* name) :
80 RooAbsPdf(other, name),
81 _x("x", this, other._x),
82 _coefList("coefList",this,other._coefList)
83{
84}
85
86////////////////////////////////////////////////////////////////////////////////
87
88/// Force use of a given normalisation range.
89/// Needed for functions or PDFs (e.g. RooAddPdf) whose shape depends on the choice of normalisation.
90void RooBernstein::selectNormalizationRange(const char* rangeName, Bool_t force)
91{
92 if (rangeName && (force || !_refRangeName.empty())) {
93 _refRangeName = rangeName;
94 }
95}
96
97////////////////////////////////////////////////////////////////////////////////
98
100{
101 double xmax,xmin;
102 std::tie(xmin, xmax) = _x->getRange(_refRangeName.empty() ? nullptr : _refRangeName.c_str());
103 Double_t x = (_x - xmin) / (xmax - xmin); // rescale to [0,1]
104 Int_t degree = _coefList.getSize() - 1; // n+1 polys of degree n
106
107 if(degree == 0) {
108
109 return ((RooAbsReal *)iter.next())->getVal();
110
111 } else if(degree == 1) {
112
113 Double_t a0 = ((RooAbsReal *)iter.next())->getVal(); // c0
114 Double_t a1 = ((RooAbsReal *)iter.next())->getVal() - a0; // c1 - c0
115 return a1 * x + a0;
116
117 } else if(degree == 2) {
118
119 Double_t a0 = ((RooAbsReal *)iter.next())->getVal(); // c0
120 Double_t a1 = 2 * (((RooAbsReal *)iter.next())->getVal() - a0); // 2 * (c1 - c0)
121 Double_t a2 = ((RooAbsReal *)iter.next())->getVal() - a1 - a0; // c0 - 2 * c1 + c2
122 return (a2 * x + a1) * x + a0;
123
124 } else if(degree > 2) {
125
126 Double_t t = x;
127 Double_t s = 1 - x;
128
129 Double_t result = ((RooAbsReal *)iter.next())->getVal() * s;
130 for(Int_t i = 1; i < degree; i++) {
131 result = (result + t * TMath::Binomial(degree, i) * ((RooAbsReal *)iter.next())->getVal()) * s;
132 t *= x;
133 }
134 result += t * ((RooAbsReal *)iter.next())->getVal();
135
136 return result;
137 }
138
139 // in case list of arguments passed is empty
140 return TMath::SignalingNaN();
141}
142
143////////////////////////////////////////////////////////////////////////////////
144/// Compute multiple values of Bernstein distribution.
145void RooBernstein::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooBatchCompute::DataMap& dataMap) const
146{
147 const int nCoef = _coefList.size();
148 std::vector<double> extraArgs(nCoef+2);
149 for (int i=0; i<nCoef; i++)
150 extraArgs[i] = static_cast<RooAbsReal&>(_coefList[i]).getVal();
151 extraArgs[nCoef] = _x.min();
152 extraArgs[nCoef+1] = _x.max();
153
155 dispatch->compute(stream, RooBatchCompute::Bernstein, output, nEvents, dataMap, {&*_x, &*_norm}, extraArgs);
156}
157
158////////////////////////////////////////////////////////////////////////////////
159
160Int_t RooBernstein::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
161{
162
163 if (matchArgs(allVars, analVars, _x)) return 1;
164 return 0;
165}
166
167////////////////////////////////////////////////////////////////////////////////
168
169Double_t RooBernstein::analyticalIntegral(Int_t code, const char* rangeName) const
170{
171 R__ASSERT(code==1) ;
172
173 double xmax,xmin;
174 std::tie(xmin, xmax) = _x->getRange(_refRangeName.empty() ? nullptr : _refRangeName.c_str());
175 const Double_t xlo = (_x.min(rangeName) - xmin) / (xmax - xmin);
176 const Double_t xhi = (_x.max(rangeName) - xmin) / (xmax - xmin);
177
178 Int_t degree= _coefList.getSize()-1; // n+1 polys of degree n
179 Double_t norm(0) ;
180
182 Double_t temp=0;
183 for (int i=0; i<=degree; ++i){
184 // for each of the i Bernstein basis polynomials
185 // represent it in the 'power basis' (the naive polynomial basis)
186 // where the integral is straight forward.
187 temp = 0;
188 for (int j=i; j<=degree; ++j){ // power basis≈ß
189 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);
190 }
191 temp *= ((RooAbsReal*)iter.next())->getVal(); // include coeff
192 norm += temp; // add this basis's contribution to total
193 }
194
195 norm *= xmax-xmin;
196 return norm;
197}
#define ClassImp(name)
Definition: Rtypes.h:364
#define R__ASSERT(e)
Definition: TError.h:118
char name[80]
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:74
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.
RooAbsReal * _norm
Definition: RooAbsPdf.h:364
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:63
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:93
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:35
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const DataMap &, const VarVector &, const ArgVector &={})=0
Bernstein basis polynomials are positive-definite in the range [0,1].
Definition: RooBernstein.h:26
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:45
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooBatchCompute::DataMap &) const
Compute multiple values of Bernstein distribution.
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
RooTemplateProxy< RooAbsRealLValue > _x
Definition: RooBernstein.h:43
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:44
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) override
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.
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
R__EXTERN RooBatchComputeInterface * dispatchCUDA
std::map< DataKey, RooSpan< const double > > DataMap
R__EXTERN RooBatchComputeInterface * dispatchCPU
This dispatch pointer points to an implementation of the compute library, provided one has been loade...
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:2085
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:685
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN)
Definition: TMath.h:858
static void output(int code)
Definition: gifencode.c:226