ROOT   Reference Guide
Searching...
No Matches
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 *
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 }
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.
188
189 RooSpan<const double> xData = _x->getValues(evalData, normSet);
190 size_t batchSize = xData.size();
191 RooSpan<double> output = evalData.makeBatch(this, batchSize);
192 const Double_t xmin = _x.min(_refRangeName?_refRangeName->GetName() : nullptr);
193 const Double_t xmax = _x.max(_refRangeName?_refRangeName->GetName() : nullptr);
194
195 const size_t nCoef = _coefList.size();
196 std::vector<double> coef(nCoef);
197 for (size_t i=0; i<nCoef; i++) {
198 coef[i] = static_cast<const RooAbsReal &>(_coefList[i]).getVal();
199 }
200 RooBatchCompute::dispatch->computeChebychev(batchSize, output.data(), xData.data(), xmin, xmax, coef);
201 return output;
202}
203
204////////////////////////////////////////////////////////////////////////////////
205
206
207Int_t RooChebychev::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /* rangeName */) const
208{
209 if (matchArgs(allVars, analVars, _x)) return 1;
210 return 0;
211}
212
213////////////////////////////////////////////////////////////////////////////////
214
215Double_t RooChebychev::analyticalIntegral(Int_t code, const char* rangeName) const
216{
217 assert(1 == code); (void)code;
218
221 const Double_t halfrange = .5 * (xmax - xmin), mid = .5 * (xmax + xmin);
222 // the full range of the function is mapped to the normalised [-1, 1] range
223 const Double_t b = (_x.max(rangeName) - mid) / halfrange;
224 const Double_t a = (_x.min(rangeName) - mid) / halfrange;
225
226 // take care to multiply with the right factor to account for the mapping to
227 // normalised range [-1, 1]
228 return halfrange * evalAnaInt(a, b);
229}
230
231////////////////////////////////////////////////////////////////////////////////
232
234{
235 // coefficient for integral(T_0(x)) is 1 (implicit), integrate by hand
236 // T_0(x) and T_1(x), and use for n > 1: integral(T_n(x) dx) =
237 // (T_n+1(x) / (n + 1) - T_n-1(x) / (n - 1)) / 2
238 double sum = b - a; // integrate T_0(x) by hand
239
240 using size_type = typename RooListProxy::Storage_t::size_type;
241 const size_type iend = _coefList.size();
242 if (iend > 0) {
243 {
244 // integrate T_1(x) by hand...
245 const double c = static_cast<const RooAbsReal &>(_coefList[0]).getVal();
246 sum = fast_fma(0.5 * (b + a) * (b - a), c, sum);
247 }
248 if (1 < iend) {
249 ChebychevIterator<double, Kind::First> bit(b), ait(a);
250 ++bit, ++ait;
251 double nminus1 = 1.;
252 for (size_type i = 1; iend != i; ++i) {
253 // integrate using recursion relation
254 const double c = static_cast<const RooAbsReal &>(_coefList[i]).getVal();
255 const double term2 = (*bit - *ait) / nminus1;
256 ++bit, ++ait, ++nminus1;
257 const double term1 = (bit.lookahead() - ait.lookahead()) / (nminus1 + 1.);
258 const double intTn = 0.5 * (term1 - term2);
259 sum = fast_fma(intTn, c, sum);
260 }
261 }
262 }
263 return sum;
264}
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define coutE(a)
double Double_t
Definition RtypesCore.h:59
#define ClassImp(name)
Definition Rtypes.h:364
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
Binding & operator=(OUT(*fun)(void))
typedef void((*Func_t)())
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:61
virtual RooSpan< const double > getValues(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet=nullptr) const
by this change, please consult the release notes for ROOT 6.24 for guidance on how to make this trans...
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:91
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:29
virtual void computeChebychev(size_t batchSize, double *__restrict output, const double *__restrict const xData, double xmin, double xmax, std::vector< double > coef)=0
Chebychev polynomial p.d.f.
RooSpan< double > evaluateSpan(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet) const
Compute multiple values of Chebychev.
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
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
TNamed * _refRangeName
RooListProxy _coefList
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.
static RooNameReg & instance()
Return reference to singleton instance.
A simple container to hold a batch of data values.
Definition RooSpan.h:34
constexpr std::span< T >::pointer data() const
Definition RooSpan.h:106
constexpr std::span< T >::index_type size() const noexcept
Definition RooSpan.h:121
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)
R__EXTERN RooBatchComputeInterface * dispatch
This dispatch pointer points to an implementation of the compute library, provided one has been loade...
This struct enables passing computation data around between elements of a computation graph.
Definition RunContext.h:31
RooSpan< double > makeBatch(const RooAbsReal *owner, std::size_t size)
Create a writable batch.
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345
static void output(int code)
Definition gifencode.c:226