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