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
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_t 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
166 // transform to range [-1, +1]
167 const Double_t 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
186namespace {
187//Author: Emmanouil Michalainas, CERN 12 AUGUST 2019
188
189void compute( size_t batchSize, double xmax, double xmin,
190 double * __restrict output,
191 const double * __restrict const xData,
192 const RooListProxy& _coefList)
193{
194 constexpr size_t block = 128;
195 const size_t nCoef = _coefList.size();
196 double prev[block][2], X[block];
197
198 for (size_t i=0; i<batchSize; i+=block) {
199 size_t stop = (i+block >= batchSize) ? batchSize-i : block;
200
201 // set a0-->prev[j][0] and a1-->prev[j][1]
202 // and x tranfsformed to range[-1..1]-->X[j]
203 for (size_t j=0; j<stop; j++) {
204 prev[j][0] = output[i+j] = 1.0;
205 prev[j][1] = X[j] = (xData[i+j] -0.5*(xmax + xmin)) / (0.5*(xmax - xmin));
206 }
207
208 for (size_t k=0; k<nCoef; k++) {
209 const double coef = static_cast<const RooAbsReal &>(_coefList[k]).getVal();
210 for (size_t j=0; j<stop; j++) {
211 output[i+j] += prev[j][1]*coef;
212
213 //compute next order
214 const double next = 2*X[j]*prev[j][1] -prev[j][0];
215 prev[j][0] = prev[j][1];
216 prev[j][1] = next;
217 }
218 }
219 }
220}
221};
222
223
224RooSpan<double> RooChebychev::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
225 auto xData = _x.getValBatch(begin, batchSize);
226 if (xData.empty()) {
227 return {};
228 }
229
230 batchSize = xData.size();
231 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
232 const Double_t xmax = _x.max(_refRangeName?_refRangeName->GetName() : nullptr);
233 const Double_t xmin = _x.min(_refRangeName?_refRangeName->GetName() : nullptr);
234 compute(batchSize, xmax, xmin, output.data(), xData.data(), _coefList);
235 return output;
236}
237////////////////////////////////////////////////////////////////////////////////
238
239Int_t RooChebychev::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /* rangeName */) const
240{
241 if (matchArgs(allVars, analVars, _x)) return 1;
242 return 0;
243}
244
245////////////////////////////////////////////////////////////////////////////////
246
247Double_t RooChebychev::analyticalIntegral(Int_t code, const char* rangeName) const
248{
249 assert(1 == code); (void)code;
250
253 const Double_t halfrange = .5 * (xmax - xmin), mid = .5 * (xmax + xmin);
254 // the full range of the function is mapped to the normalised [-1, 1] range
255 const Double_t b = (_x.max(rangeName) - mid) / halfrange;
256 const Double_t a = (_x.min(rangeName) - mid) / halfrange;
257
258 // take care to multiply with the right factor to account for the mapping to
259 // normalised range [-1, 1]
260 return halfrange * evalAnaInt(a, b);
261}
262
263////////////////////////////////////////////////////////////////////////////////
264
266{
267 // coefficient for integral(T_0(x)) is 1 (implicit), integrate by hand
268 // T_0(x) and T_1(x), and use for n > 1: integral(T_n(x) dx) =
269 // (T_n+1(x) / (n + 1) - T_n-1(x) / (n - 1)) / 2
270 double sum = b - a; // integrate T_0(x) by hand
271
272 using size_type = typename RooListProxy::Storage_t::size_type;
273 const size_type iend = _coefList.size();
274 if (iend > 0) {
275 {
276 // integrate T_1(x) by hand...
277 const double c = static_cast<const RooAbsReal &>(_coefList[0]).getVal();
278 sum = fast_fma(0.5 * (b + a) * (b - a), c, sum);
279 }
280 if (1 < iend) {
281 ChebychevIterator<double, Kind::First> bit(b), ait(a);
282 ++bit, ++ait;
283 double nminus1 = 1.;
284 for (size_type i = 1; iend != i; ++i) {
285 // integrate using recursion relation
286 const double c = static_cast<const RooAbsReal &>(_coefList[i]).getVal();
287 const double term2 = (*bit - *ait) / nminus1;
288 ++bit, ++ait, ++nminus1;
289 const double term1 = (bit.lookahead() - ait.lookahead()) / (nminus1 + 1.);
290 const double intTn = 0.5 * (term1 - term2);
291 sum = fast_fma(intTn, c, sum);
292 }
293 }
294 }
295 return sum;
296}
#define b(i)
Definition: RSha256.hxx:100
#define c(i)
Definition: RSha256.hxx:101
#define coutE(a)
Definition: RooMsgService.h:34
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
Binding & operator=(OUT(*fun)(void))
typedef void((*Func_t)())
TTime operator*(const TTime &t1, const TTime &t2)
Definition: TTime.h:85
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:58
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:59
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:87
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:447
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Chebychev polynomial p.d.f.
Definition: RooChebychev.h:25
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
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:44
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
TNamed * _refRangeName
Definition: RooChebychev.h:46
RooListProxy _coefList
Definition: RooChebychev.h:45
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.
RooListProxy is the concrete proxy for RooArgList objects.
Definition: RooListProxy.h:25
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
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
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
Double_t min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
RooSpan< const double > getValBatch(std::size_t begin, std::size_t batchSize) const
Double_t 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
constexpr size_t block
Definition: BatchHelpers.h:29
double T(double x)
Definition: ChebyshevPol.h:34
@ InputArguments
Definition: RooGlobalFunc.h:68
auto * a
Definition: textangle.C:12
static long int sum(long int i)
Definition: Factory.cxx:2276
static void output(int code)
Definition: gifencode.c:226