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