Logo ROOT  
Reference Guide
RooChebychev.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * GR, Gerhard Raven, UC San Diego, Gerhard.Raven@slac.stanford.edu
7 * *
8 * Copyright (c) 2000-2005, Regents of the University of California *
9 * and Stanford University. All rights reserved. *
10 * *
11 * Redistribution and use in source and binary forms, *
12 * with or without modification, are permitted according to the terms *
13 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
14 *****************************************************************************/
15
16/** \class RooChebychev
17 \ingroup Roofit
18
19Chebychev polynomial p.d.f. of the first kind.
20
21The coefficient that goes with \f$ T_0(x)=1 \f$ (i.e. the constant polynomial) is
22implicitly assumed to be 1, and the list of coefficients supplied by callers
23starts with the coefficient that goes with \f$ T_1(x)=x \f$ (i.e. the linear term).
24**/
25
26#include "RooChebychev.h"
27#include "RooFit.h"
28#include "RooAbsReal.h"
29#include "RooRealVar.h"
30#include "RooArgList.h"
31#include "RooNameReg.h"
32#include "RooBatchCompute.h"
33
34#include <cmath>
35
37
38namespace { // anonymous namespace to hide implementation details
39/// use fast FMA if available, fall back to normal arithmetic if not
40static inline double fast_fma(
41 const double x, const double y, const double z) noexcept
42{
43#if defined(FP_FAST_FMA) // check if std::fma has fast hardware implementation
44 return std::fma(x, y, z);
45#else // defined(FP_FAST_FMA)
46 // std::fma might be slow, so use a more pedestrian implementation
47#if defined(__clang__)
48#pragma STDC FP_CONTRACT ON // hint clang that using an FMA is okay here
49#endif // defined(__clang__)
50 return (x * y) + z;
51#endif // defined(FP_FAST_FMA)
52}
53
54/// Chebychev polynomials of first or second kind
55enum class Kind : int { First = 1, Second = 2 };
56
57/** @brief ChebychevIterator evaluates increasing orders at given x
58 *
59 * @author Manuel Schiller <Manuel.Schiller@glasgow.ac.uk>
60 * @date 2019-03-24
61 */
62template <typename T, Kind KIND>
63class ChebychevIterator {
64private:
65 T _last = 1;
66 T _curr = 0;
67 T _twox = 0;
68
69public:
70 /// default constructor
71 constexpr ChebychevIterator() = default;
72 /// copy constructor
73 ChebychevIterator(const ChebychevIterator &) = default;
74 /// move constructor
75 ChebychevIterator(ChebychevIterator &&) = default;
76 /// construct from given x in [-1, 1]
77 constexpr ChebychevIterator(const T &x)
78 : _curr(static_cast<int>(KIND) * x), _twox(2 * x)
79 {}
80
81 /// (copy) assignment
82 ChebychevIterator &operator=(const ChebychevIterator &) = default;
83 /// move assignment
84 ChebychevIterator &operator=(ChebychevIterator &&) = default;
85
86 /// get value of Chebychev polynomial at current order
87 constexpr inline T operator*() const noexcept { return _last; }
88 // get value of Chebychev polynomial at (current + 1) order
89 constexpr inline T lookahead() const noexcept { return _curr; }
90 /// move on to next order, return reference to new value
91 inline ChebychevIterator &operator++() noexcept
92 {
93 //T newval = fast_fma(_twox, _curr, -_last);
94 T newval = _twox*_curr -_last;
95 _last = _curr;
96 _curr = newval;
97 return *this;
98 }
99 /// move on to next order, return copy of new value
100 inline ChebychevIterator operator++(int) noexcept
101 {
102 ChebychevIterator retVal(*this);
103 operator++();
104 return retVal;
105 }
106};
107} // anonymous namespace
108
109////////////////////////////////////////////////////////////////////////////////
110
111RooChebychev::RooChebychev() : _refRangeName(0)
112{
113}
114
115////////////////////////////////////////////////////////////////////////////////
116/// Constructor
117
118RooChebychev::RooChebychev(const char* name, const char* title,
119 RooAbsReal& x, const RooArgList& coefList):
120 RooAbsPdf(name, title),
121 _x("x", "Dependent", this, x),
122 _coefList("coefficients","List of coefficients",this),
123 _refRangeName(0)
124{
125 for (const auto coef : coefList) {
126 if (!dynamic_cast<RooAbsReal*>(coef)) {
127 coutE(InputArguments) << "RooChebychev::ctor(" << GetName() <<
128 ") ERROR: coefficient " << coef->GetName() <<
129 " is not of type RooAbsReal" << std::endl ;
130 throw std::invalid_argument("Wrong input arguments for RooChebychev");
131 }
132 _coefList.add(*coef) ;
133 }
134}
135
136////////////////////////////////////////////////////////////////////////////////
137
138RooChebychev::RooChebychev(const RooChebychev& other, const char* name) :
139 RooAbsPdf(other, name),
140 _x("x", this, other._x),
141 _coefList("coefList",this,other._coefList),
142 _refRangeName(other._refRangeName)
143{
144}
145
146////////////////////////////////////////////////////////////////////////////////
147
148void RooChebychev::selectNormalizationRange(const char* rangeName, Bool_t force)
149{
150 if (rangeName && (force || !_refRangeName)) {
152 }
153 if (!rangeName) {
154 _refRangeName = 0 ;
155 }
156}
157
158////////////////////////////////////////////////////////////////////////////////
159
161{
162 // first bring the range of the variable _x to the normalised range [-1, 1]
163 // calculate sum_k c_k T_k(x) where x is given in the normalised range,
164 // c_0 = 1, and the higher coefficients are given in _coefList
167 // transform to range [-1, +1]
168 const Double_t x = (_x - 0.5 * (xmax + xmin)) / (0.5 * (xmax - xmin));
169 // extract current values of coefficients
170 using size_type = typename RooListProxy::Storage_t::size_type;
171 const size_type iend = _coefList.size();
172 double sum = 1.;
173 if (iend > 0) {
174 ChebychevIterator<double, Kind::First> chit(x);
175 ++chit;
176 for (size_type i = 0; iend != i; ++i, ++chit) {
177 auto c = static_cast<const RooAbsReal &>(_coefList[i]).getVal();
178 //sum = fast_fma(*chit, c, sum);
179 sum += *chit*c;
180 }
181 }
182 return sum;
183}
184
185////////////////////////////////////////////////////////////////////////////////
186/// Compute multiple values of Chebychev.
187void RooChebychev::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooBatchCompute::DataMap& dataMap) const
188{
190 for (auto* coef:_coefList)
191 extraArgs.push_back( static_cast<const RooAbsReal*>(coef)->getVal() );
192 extraArgs.push_back( _x.min(_refRangeName?_refRangeName->GetName() : nullptr) );
193 extraArgs.push_back( _x.max(_refRangeName?_refRangeName->GetName() : nullptr) );
195 dispatch->compute(stream, RooBatchCompute::Chebychev, output, nEvents, dataMap, {&*_x,&*_norm}, extraArgs);
196}
197
198////////////////////////////////////////////////////////////////////////////////
199
200
201Int_t RooChebychev::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /* rangeName */) const
202{
203 if (matchArgs(allVars, analVars, _x)) return 1;
204 return 0;
205}
206
207////////////////////////////////////////////////////////////////////////////////
208
209Double_t RooChebychev::analyticalIntegral(Int_t code, const char* rangeName) const
210{
211 assert(1 == code); (void)code;
212
215 const Double_t halfrange = .5 * (xmax - xmin), mid = .5 * (xmax + xmin);
216 // the full range of the function is mapped to the normalised [-1, 1] range
217 const Double_t b = (_x.max(rangeName) - mid) / halfrange;
218 const Double_t a = (_x.min(rangeName) - mid) / halfrange;
219
220 // take care to multiply with the right factor to account for the mapping to
221 // normalised range [-1, 1]
222 return halfrange * evalAnaInt(a, b);
223}
224
225////////////////////////////////////////////////////////////////////////////////
226
228{
229 // coefficient for integral(T_0(x)) is 1 (implicit), integrate by hand
230 // T_0(x) and T_1(x), and use for n > 1: integral(T_n(x) dx) =
231 // (T_n+1(x) / (n + 1) - T_n-1(x) / (n - 1)) / 2
232 double sum = b - a; // integrate T_0(x) by hand
233
234 using size_type = typename RooListProxy::Storage_t::size_type;
235 const size_type iend = _coefList.size();
236 if (iend > 0) {
237 {
238 // integrate T_1(x) by hand...
239 const double c = static_cast<const RooAbsReal &>(_coefList[0]).getVal();
240 sum = fast_fma(0.5 * (b + a) * (b - a), c, sum);
241 }
242 if (1 < iend) {
243 ChebychevIterator<double, Kind::First> bit(b), ait(a);
244 ++bit, ++ait;
245 double nminus1 = 1.;
246 for (size_type i = 1; iend != i; ++i) {
247 // integrate using recursion relation
248 const double c = static_cast<const RooAbsReal &>(_coefList[i]).getVal();
249 const double term2 = (*bit - *ait) / nminus1;
250 ++bit, ++ait, ++nminus1;
251 const double term1 = (bit.lookahead() - ait.lookahead()) / (nminus1 + 1.);
252 const double intTn = 0.5 * (term1 - term2);
253 sum = fast_fma(intTn, c, sum);
254 }
255 }
256 }
257 return sum;
258}
typedef void(GLAPIENTRYP _GLUfuncptr)(void)
#define b(i)
Definition: RSha256.hxx:100
#define c(i)
Definition: RSha256.hxx:101
#define coutE(a)
Definition: RooMsgService.h:33
double Double_t
Definition: RtypesCore.h:59
#define ClassImp(name)
Definition: Rtypes.h:364
char name[80]
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
Binding & operator=(OUT(*fun)(void))
TTime operator*(const TTime &t1, const TTime &t2)
Definition: TTime.h:85
Storage_t::size_type size() const
RooAbsReal * _norm
Definition: RooAbsPdf.h:364
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
Chebychev polynomial p.d.f.
Definition: RooChebychev.h:25
virtual void selectNormalizationRange(const char *rangeName=0, Bool_t force=kFALSE)
Interface function to force use of a given normalization range to interpret function value.
RooRealProxy _x
Definition: RooChebychev.h:43
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooBatchCompute::DataMap &) const
Compute multiple values of Chebychev.
TNamed * _refRangeName
Definition: RooChebychev.h:45
RooListProxy _coefList
Definition: RooChebychev.h:44
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Double_t evalAnaInt(const Double_t a, const Double_t b) const
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Reimplementation of standard RooArgList::add()
const TNamed * constPtr(const char *stringPtr)
Return a unique TNamed pointer for given C++ string.
Definition: RooNameReg.cxx:61
static RooNameReg & instance()
Return reference to singleton instance.
Definition: RooNameReg.cxx:51
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.
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
double T(double x)
Definition: ChebyshevPol.h:34
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...
std::vector< double > ArgVector
@ InputArguments
Definition: RooGlobalFunc.h:61
auto * a
Definition: textangle.C:12
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345
static void output(int code)
Definition: gifencode.c:226