Logo ROOT   master
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 
19 Chebychev polynomial p.d.f. of the first kind.
20 
21 The coefficient that goes with \f$ T_0(x)=1 \f$ (i.e. the constant polynomial) is
22 implicitly assumed to be 1, and the list of coefficients supplied by callers
23 starts 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 
37 namespace { // anonymous namespace to hide implementation details
38 /// use fast FMA if available, fall back to normal arithmetic if not
39 static 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
54 enum 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  */
61 template <typename T, Kind KIND>
62 class ChebychevIterator {
63 private:
64  T _last = 1;
65  T _curr = 0;
66  T _twox = 0;
67 
68 public:
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 
110 RooChebychev::RooChebychev() : _refRangeName(0)
111 {
112 }
113 
114 ////////////////////////////////////////////////////////////////////////////////
115 /// Constructor
116 
117 RooChebychev::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 
137 RooChebychev::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 
147 void 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 
186 namespace {
187 //Author: Emmanouil Michalainas, CERN 12 AUGUST 2019
188 
189 void 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 
224 RooSpan<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 
239 Int_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 
247 Double_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 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
double max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
static long int sum(long int i)
Definition: Factory.cxx:2272
#define coutE(a)
Definition: RooMsgService.h:33
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
float xmin
Definition: THbookFile.cxx:93
Storage_t::size_type size() const
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t mid
Definition: TRolke.cxx:630
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
double T(double x)
Definition: ChebyshevPol.h:34
Kind
Chebychev polynomials of first or second kind.
RooRealProxy _x
Definition: RooChebychev.h:44
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral. ...
double min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
Double_t evalAnaInt(const Double_t a, const Double_t b) const
Double_t x[n]
Definition: legend1.C:17
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
RooSpan< const typename T::value_type > getValBatch(std::size_t begin, std::size_t batchSize) const
Retrieve a batch of real or category data.
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize, const RooArgSet *const normSet=nullptr, Tag_t ownerTag=kUnspecified)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:118
TNamed * _refRangeName
Definition: RooChebychev.h:46
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
RooListProxy _coefList
Definition: RooChebychev.h:45
static RooNameReg & instance()
Return reference to singleton instance.
Definition: RooNameReg.cxx:51
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
auto * a
Definition: textangle.C:12
TTime operator*(const TTime &t1, const TTime &t2)
Definition: TTime.h:85
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported...
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
float xmax
Definition: THbookFile.cxx:93
RooListProxy is the concrete proxy for RooArgList objects.
Definition: RooListProxy.h:25
#define ClassImp(name)
Definition: Rtypes.h:361
double Double_t
Definition: RtypesCore.h:57
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
Double_t y[n]
Definition: legend1.C:17
Binding & operator=(OUT(*fun)(void))
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
typedef void((*Func_t)())
RooAbsPdf, the base class of all PDFs
Definition: RooAbsPdf.h:40
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...
const TNamed * constPtr(const char *stringPtr)
Return a unique TNamed pointer for given C++ string.
Definition: RooNameReg.cxx:61
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
#define c(i)
Definition: RSha256.hxx:101
constexpr size_t block
Definition: BatchHelpers.h:29
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
Chebychev polynomial p.d.f.
Definition: RooChebychev.h:25
static void output(int code)
Definition: gifencode.c:226
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:450
char name[80]
Definition: TGX11.cxx:109