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 "RooAbsReal.h"
37#include "RooRealVar.h"
38#include "RooArgList.h"
39#include "RooBatchCompute.h"
40
41#include "TMath.h"
42
43#include <cmath>
44using namespace std;
45
47
48/// Constructor
49////////////////////////////////////////////////////////////////////////////////
50
51RooBernstein::RooBernstein(const char* name, const char* title,
52 RooAbsRealLValue& x, const RooArgList& coefList):
53 RooAbsPdf(name, title),
54 _x("x", "Dependent", this, x),
55 _coefList("coefficients","List of coefficients",this)
56{
57 for (auto *coef : coefList) {
58 if (!dynamic_cast<RooAbsReal*>(coef)) {
59 cout << "RooBernstein::ctor(" << GetName() << ") ERROR: coefficient " << coef->GetName()
60 << " is not of type RooAbsReal" << endl ;
61 R__ASSERT(0) ;
62 }
64 }
65}
66
67////////////////////////////////////////////////////////////////////////////////
68
70 : RooAbsPdf(other, name), _x("x", this, other._x), _coefList("coefList", this, other._coefList)
71{
72}
73
74////////////////////////////////////////////////////////////////////////////////
75
76/// Force use of a given normalisation range.
77/// Needed for functions or PDFs (e.g. RooAddPdf) whose shape depends on the choice of normalisation.
78void RooBernstein::selectNormalizationRange(const char* rangeName, bool force)
79{
80 if (rangeName && (force || !_refRangeName.empty())) {
81 _refRangeName = rangeName;
82 }
83}
84
85////////////////////////////////////////////////////////////////////////////////
86
88{
89 double xmax,xmin;
90 std::tie(xmin, xmax) = _x->getRange(_refRangeName.empty() ? nullptr : _refRangeName.c_str());
91 double x = (_x - xmin) / (xmax - xmin); // rescale to [0,1]
92 Int_t degree = _coefList.getSize() - 1; // n+1 polys of degree n
93
94 if(degree == 0) {
95
96 return static_cast<RooAbsReal &>(_coefList[0]).getVal();
97
98 } else if(degree == 1) {
99
100 double a0 = static_cast<RooAbsReal &>(_coefList[0]).getVal(); // c0
101 double a1 = static_cast<RooAbsReal &>(_coefList[1]).getVal() - a0; // c1 - c0
102 return a1 * x + a0;
103
104 } else if(degree == 2) {
105
106 double a0 = static_cast<RooAbsReal &>(_coefList[0]).getVal(); // c0
107 double a1 = 2 * (static_cast<RooAbsReal &>(_coefList[1]).getVal() - a0); // 2 * (c1 - c0)
108 double a2 = static_cast<RooAbsReal &>(_coefList[2]).getVal() - a1 - a0; // c0 - 2 * c1 + c2
109 return (a2 * x + a1) * x + a0;
110
111 } else if(degree > 2) {
112
113 double t = x;
114 double s = 1 - x;
115
116 double result = static_cast<RooAbsReal &>(_coefList[0]).getVal() * s;
117 for(Int_t i = 1; i < degree; i++) {
118 result = (result + t * TMath::Binomial(degree, i)
119 * static_cast<RooAbsReal &>(_coefList[i]).getVal()) * s;
120 t *= x;
121 }
122 result += t * static_cast<RooAbsReal &>(_coefList[degree]).getVal();
123
124 return result;
125 }
126
127 // in case list of arguments passed is empty
128 return TMath::SignalingNaN();
129}
130
131////////////////////////////////////////////////////////////////////////////////
132/// Compute multiple values of Bernstein distribution.
133void RooBernstein::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
134{
135 const int nCoef = _coefList.size();
136 std::vector<double> extraArgs(nCoef+2);
137 for (int i=0; i<nCoef; i++)
138 extraArgs[i] = static_cast<RooAbsReal&>(_coefList[i]).getVal();
139 extraArgs[nCoef] = _x.min();
140 extraArgs[nCoef+1] = _x.max();
141
143 dispatch->compute(stream, RooBatchCompute::Bernstein, output, nEvents, {dataMap.at(_x)}, extraArgs);
144}
145
146////////////////////////////////////////////////////////////////////////////////
147
148Int_t RooBernstein::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
149{
150
151 if (matchArgs(allVars, analVars, _x)) return 1;
152 return 0;
153}
154
155////////////////////////////////////////////////////////////////////////////////
156
157double RooBernstein::analyticalIntegral(Int_t code, const char* rangeName) const
158{
159 R__ASSERT(code==1) ;
160
161 double xmax,xmin;
162 std::tie(xmin, xmax) = _x->getRange(_refRangeName.empty() ? nullptr : _refRangeName.c_str());
163 const double xlo = (_x.min(rangeName) - xmin) / (xmax - xmin);
164 const double xhi = (_x.max(rangeName) - xmin) / (xmax - xmin);
165
166 Int_t degree= _coefList.getSize()-1; // n+1 polys of degree n
167 double norm(0) ;
168
169 double temp=0;
170 for (int i=0; i<=degree; ++i){
171 // for each of the i Bernstein basis polynomials
172 // represent it in the 'power basis' (the naive polynomial basis)
173 // where the integral is straight forward.
174 temp = 0;
175 for (int j=i; j<=degree; ++j){ // power basis≈ß
176 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);
177 }
178 temp *= static_cast<RooAbsReal &>(_coefList[i]).getVal(); // include coeff
179 norm += temp; // add this basis's contribution to total
180 }
181
182 norm *= xmax-xmin;
183 return norm;
184}
#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:64
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:94
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:57
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
std::string _refRangeName
Definition: RooBernstein.h:46
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
double analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
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.
void selectNormalizationRange(const char *rangeName=0, bool force=false) override
Force use of a given normalisation range.
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 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.
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:1753
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:2092
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:718
Double_t SignalingNaN()
Returns a signaling NaN as defined by IEEE 754](http://en.wikipedia.org/wiki/NaN#Signaling_NaN).
Definition: TMath.h:907
static void output()