ROOT   Reference Guide
Searching...
No Matches
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.
146 RooSpan<const double> xData = _x->getValues(evalData, normSet);
147 const size_t batchSize = xData.size();
148 RooSpan<double> output = evalData.makeBatch(this, batchSize);
149
150 const int nCoef = _coefList.size();
151 std::vector<double> coef(nCoef);
152 for (int i=0; i<nCoef; i++) {
153 coef[i] = static_cast<RooAbsReal&>(_coefList[i]).getVal();
154 }
155 RooBatchCompute::dispatch->computeBernstein(batchSize, output.data(), xData.data(), _x.min(), _x.max(), coef);
156 return output;
157}
158
159////////////////////////////////////////////////////////////////////////////////
160
161Int_t RooBernstein::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
162{
163
164 if (matchArgs(allVars, analVars, _x)) return 1;
165 return 0;
166}
167
168////////////////////////////////////////////////////////////////////////////////
169
170Double_t RooBernstein::analyticalIntegral(Int_t code, const char* rangeName) const
171{
172 R__ASSERT(code==1) ;
173
174 double xmax,xmin;
175 std::tie(xmin, xmax) = _x->getRange(_refRangeName.empty() ? nullptr : _refRangeName.c_str());
176 const Double_t xlo = (_x.min(rangeName) - xmin) / (xmax - xmin);
177 const Double_t xhi = (_x.max(rangeName) - xmin) / (xmax - xmin);
178
179 Int_t degree= _coefList.getSize()-1; // n+1 polys of degree n
180 Double_t norm(0) ;
181
183 Double_t temp=0;
184 for (int i=0; i<=degree; ++i){
185 // for each of the i Bernstein basis polynomials
186 // represent it in the 'power basis' (the naive polynomial basis)
187 // where the integral is straight forward.
188 temp = 0;
189 for (int j=i; j<=degree; ++j){ // power basis≈ß
190 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);
191 }
192 temp *= ((RooAbsReal*)iter.next())->getVal(); // include coeff
193 norm += temp; // add this basis's contribution to total
194 }
195
196 norm *= xmax-xmin;
197 return norm;
198}
#define ClassImp(name)
Definition Rtypes.h:364
#define R__ASSERT(e)
Definition TError.h:120
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
double pow(double, double)
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition RooAbsArg.h:72
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:61
virtual RooSpan< const double > getValues(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet=nullptr) const
by this change, please consult the release notes for ROOT 6.24 for guidance on how to make this trans...
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:91
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:29
virtual void computeBernstein(size_t batchSize, double *__restrict output, const double *__restrict const xData, double xmin, double xmax, std::vector< double > coef)=0
Bernstein basis polynomials are positive-definite in the range [0,1].
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
std::string _refRangeName
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
RooSpan< double > evaluateSpan(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet) const
Compute multiple values of Bernstein distribution.
RooTemplateProxy< RooAbsRealLValue > _x
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
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
A simple container to hold a batch of data values.
Definition RooSpan.h:34
constexpr std::span< T >::pointer data() const
Definition RooSpan.h:106
constexpr std::span< T >::index_type size() const noexcept
Definition RooSpan.h:121
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 * dispatch
This dispatch pointer points to an implementation of the compute library, provided one has been loade...
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:735
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN)
Definition TMath.h:908
This struct enables passing computation data around between elements of a computation graph.
Definition RunContext.h:31
RooSpan< double > makeBatch(const RooAbsReal *owner, std::size_t size)
Create a writable batch.
static void output(int code)
Definition gifencode.c:226