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 "RooRealVar.h"
37#include "RooArgList.h"
38#include "RooBatchCompute.h"
39
40#include "TMath.h"
41
42#include <cmath>
43using namespace std;
44
46
47/// Constructor
48////////////////////////////////////////////////////////////////////////////////
49
50RooBernstein::RooBernstein(const char* name, const char* title,
51 RooAbsRealLValue& x, const RooArgList& coefList):
52 RooAbsPdf(name, title),
53 _x("x", "Dependent", this, x),
54 _coefList("coefficients","List of coefficients",this)
55{
56 for (auto *coef : coefList) {
57 if (!dynamic_cast<RooAbsReal*>(coef)) {
58 cout << "RooBernstein::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName()
59 << " is not of type RooAbsReal" << endl ;
60 R__ASSERT(0) ;
61 }
62 _coefList.add(*coef) ;
63 }
64}
65
66////////////////////////////////////////////////////////////////////////////////
67
69 : RooAbsPdf(other, name), _x("x", this, other._x), _coefList("coefList", this, other._coefList)
70{
71}
72
73////////////////////////////////////////////////////////////////////////////////
74
75/// Force use of a given normalisation range.
76/// Needed for functions or PDFs (e.g. RooAddPdf) whose shape depends on the choice of normalisation.
77void RooBernstein::selectNormalizationRange(const char* rangeName, bool force)
78{
79 if (rangeName && (force || !_refRangeName.empty())) {
80 _refRangeName = rangeName;
81 }
82}
83
84////////////////////////////////////////////////////////////////////////////////
85
87{
88 double xmax,xmin;
89 std::tie(xmin, xmax) = _x->getRange(_refRangeName.empty() ? nullptr : _refRangeName.c_str());
90 double x = (_x - xmin) / (xmax - xmin); // rescale to [0,1]
91 Int_t degree = _coefList.getSize() - 1; // n+1 polys of degree n
92
93 if(degree == 0) {
94
95 return static_cast<RooAbsReal &>(_coefList[0]).getVal();
96
97 } else if(degree == 1) {
98
99 double a0 = static_cast<RooAbsReal &>(_coefList[0]).getVal(); // c0
100 double a1 = static_cast<RooAbsReal &>(_coefList[1]).getVal() - a0; // c1 - c0
101 return a1 * x + a0;
102
103 } else if(degree == 2) {
104
105 double a0 = static_cast<RooAbsReal &>(_coefList[0]).getVal(); // c0
106 double a1 = 2 * (static_cast<RooAbsReal &>(_coefList[1]).getVal() - a0); // 2 * (c1 - c0)
107 double a2 = static_cast<RooAbsReal &>(_coefList[2]).getVal() - a1 - a0; // c0 - 2 * c1 + c2
108 return (a2 * x + a1) * x + a0;
109
110 } else if(degree > 2) {
111
112 double t = x;
113 double s = 1 - x;
114
115 double result = static_cast<RooAbsReal &>(_coefList[0]).getVal() * s;
116 for(Int_t i = 1; i < degree; i++) {
117 result = (result + t * TMath::Binomial(degree, i)
118 * static_cast<RooAbsReal &>(_coefList[i]).getVal()) * s;
119 t *= x;
120 }
121 result += t * static_cast<RooAbsReal &>(_coefList[degree]).getVal();
122
123 return result;
124 }
125
126 // in case list of arguments passed is empty
127 return TMath::SignalingNaN();
128}
129
130////////////////////////////////////////////////////////////////////////////////
131/// Compute multiple values of Bernstein distribution.
132void RooBernstein::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
133{
134 const int nCoef = _coefList.size();
135 std::vector<double> extraArgs(nCoef+2);
136 for (int i=0; i<nCoef; i++)
137 extraArgs[i] = static_cast<RooAbsReal&>(_coefList[i]).getVal();
138 extraArgs[nCoef] = _x.min();
139 extraArgs[nCoef+1] = _x.max();
140
142 dispatch->compute(stream, RooBatchCompute::Bernstein, output, nEvents, {dataMap.at(_x)}, extraArgs);
143}
144
145////////////////////////////////////////////////////////////////////////////////
146
147Int_t RooBernstein::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
148{
149
150 if (matchArgs(allVars, analVars, _x)) return 1;
151 return 0;
152}
153
154////////////////////////////////////////////////////////////////////////////////
155
156double RooBernstein::analyticalIntegral(Int_t code, const char* rangeName) const
157{
158 R__ASSERT(code==1) ;
159
160 double xmax,xmin;
161 std::tie(xmin, xmax) = _x->getRange(_refRangeName.empty() ? nullptr : _refRangeName.c_str());
162 const double xlo = (_x.min(rangeName) - xmin) / (xmax - xmin);
163 const double xhi = (_x.max(rangeName) - xmin) / (xmax - xmin);
164
165 Int_t degree= _coefList.getSize()-1; // n+1 polys of degree n
166 double norm(0) ;
167
168 double temp=0;
169 for (int i=0; i<=degree; ++i){
170 // for each of the i Bernstein basis polynomials
171 // represent it in the 'power basis' (the naive polynomial basis)
172 // where the integral is straight forward.
173 temp = 0;
174 for (int j=i; j<=degree; ++j){ // power basisŧ
175 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);
176 }
177 temp *= static_cast<RooAbsReal &>(_coefList[i]).getVal(); // include coeff
178 norm += temp; // add this basis's contribution to total
179 }
180
181 norm *= xmax-xmin;
182 return norm;
183}
#define ClassImp(name)
Definition: Rtypes.h:375
#define R__ASSERT(e)
Definition: TError.h:118
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
char name[80]
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
Int_t getSize() const
Return the number of elements in the collection.
Storage_t::size_type size() const
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:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
bool matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
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:56
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const VarVector &, const ArgVector &={})=0
Bernstein basis polynomials are positive-definite in the range [0,1].
Definition: RooBernstein.h:26
void selectNormalizationRange(const char *rangeName=nullptr, bool force=false) override
Force use of a given normalisation range.
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
std::string _refRangeName
Definition: RooBernstein.h:46
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooTemplateProxy< RooAbsRealLValue > _x
Definition: RooBernstein.h:44
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Compute multiple values of Bernstein distribution.
RooListProxy _coefList
Definition: RooBernstein.h:45
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
auto & at(RooAbsArg const *arg, RooAbsArg const *=nullptr)
Definition: DataMap.h:88
double max(const char *rname=nullptr) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
double min(const char *rname=nullptr) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1770
Double_t x[n]
Definition: legend1.C:17
R__EXTERN RooBatchComputeInterface * dispatchCUDA
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)
Calculates the binomial coefficient n over k.
Definition: TMath.cxx:2110
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:719
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN).
Definition: TMath.h:908
static void output()