Skip to content

Commit 2525d2d

Browse files
fweigfweig
committed
TPCFastTransformation: Resolve recursion at compile time with templates.
1 parent bf8eb6f commit 2525d2d

File tree

4 files changed

+115
-45
lines changed

4 files changed

+115
-45
lines changed

‎GPU/Common/GPUCommonDefAPI.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -79,7 +79,7 @@
7979
#define GPUdDefault()
8080
#define GPUhdDefault()
8181
#define GPUdi() inline
82-
#define GPUdii() inline
82+
#define GPUdii() __attribute__((always_inline)) inline
8383
#define GPUdni()
8484
#define GPUdnii()
8585
#define GPUh() INVALID_TRIGGER_ERROR_NO_HOST_CODE

‎GPU/GPUTracking/DataTypes/CalibdEdxTrackTopologyPol.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,10 @@ class CalibdEdxTrackTopologyPol : public o2::gpu::FlatObject
6262
/// \param region region of the TPC
6363
/// \param charge correction for maximum or total charge
6464
/// \param x coordinates where the correction is evaluated
65-
GPUd() float getCorrection(const int32_t region, const ChargeType charge, float x[/*inpXdim*/]) const { return (charge == ChargeType::Tot) ? mCalibPolsqTot[region].eval(x) : mCalibPolsqMax[region].eval(x); }
65+
GPUd() float getCorrection(const int32_t region, const ChargeType charge, float x[/*inpXdim*/]) const
66+
{
67+
return (charge == ChargeType::Tot) ? mCalibPolsqTot[region].eval(x) : mCalibPolsqMax[region].eval(x);
68+
}
6669

6770
/// \return returns the track topology correction
6871
/// \param region region of the TPC

‎GPU/TPCFastTransformation/MultivariatePolynomialHelper.h‎

Lines changed: 92 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -60,28 +60,57 @@ class MultivariatePolynomialParametersHelper
6060
/// \returns number of parameters for given dimension and degree of polynomials
6161
/// calculates the number of parameters for a multivariate polynomial for given degree: nParameters = (n+d-1 d) -> binomial coefficient
6262
/// see: https://mathoverflow.net/questions/225953/number-of-polynomial-terms-for-certain-degree-and-certain-number-of-variables
63-
GPUd() static constexpr uint32_t getNParametersAllTerms(const uint32_t degree, const uint32_t dim) { return (degree == 0) ? binomialCoeff(dim - 1, 0) : binomialCoeff(dim - 1 + degree, degree) + getNParametersAllTerms(degree - 1, dim); }
63+
template <uint32_t Degree, uint32_t Dim>
64+
GPUd() static constexpr uint32_t getNParametersAllTerms()
65+
{
66+
if constexpr (Degree == 0) {
67+
return binomialCoeff<Dim - 1, 0>();
68+
} else {
69+
return binomialCoeff<Dim - 1 + Degree, Degree>() + getNParametersAllTerms<Degree - 1, Dim>();
70+
}
71+
}
6472

6573
/// \returns the number of parameters for interaction terms only (see: https://en.wikipedia.org/wiki/Combination)
66-
GPUd() static constexpr uint32_t getNParametersInteractionOnly(const uint32_t degree, const uint32_t dim) { return (degree == 0) ? binomialCoeff(dim - 1, 0) : binomialCoeff(dim, degree) + getNParametersInteractionOnly(degree - 1, dim); }
74+
template <uint32_t Degree, uint32_t Dim>
75+
GPUd() static constexpr uint32_t getNParametersInteractionOnly()
76+
{
77+
if constexpr (Degree == 0) {
78+
return binomialCoeff<Dim - 1, 0>();
79+
} else {
80+
return binomialCoeff<Dim, Degree>() + getNParametersInteractionOnly<Degree - 1, Dim>();
81+
}
82+
}
6783

68-
GPUd() static constexpr uint32_t getNParameters(const uint32_t degree, const uint32_t dim, const bool interactionOnly)
84+
template <uint32_t Degree, uint32_t Dim, bool InteractionOnly>
85+
GPUd() static constexpr uint32_t getNParameters()
6986
{
70-
if (interactionOnly) {
71-
return getNParametersInteractionOnly(degree, dim);
87+
if constexpr (InteractionOnly) {
88+
return getNParametersInteractionOnly<Degree, Dim>();
7289
} else {
73-
return getNParametersAllTerms(degree, dim);
90+
return getNParametersAllTerms<Degree, Dim>();
7491
}
7592
}
7693

7794
private:
7895
/// calculate factorial of n
7996
/// \return returns n!
80-
GPUd() static constexpr uint32_t factorial(const uint32_t n) { return (n == 0) || (n == 1) ? 1 : n * factorial(n - 1); }
97+
template <uint32_t N>
98+
GPUd() static constexpr uint32_t factorial()
99+
{
100+
if constexpr (N == 0 || N == 1) {
101+
return 1;
102+
} else {
103+
return N * factorial<N - 1>();
104+
}
105+
}
81106

82107
/// calculates binomial coefficient
83108
/// \return returns (n k)
84-
GPUd() static constexpr uint32_t binomialCoeff(const uint32_t n, const uint32_t k) { return factorial(n) / (factorial(k) * factorial(n - k)); }
109+
template <uint32_t N, uint32_t K>
110+
GPUd() static constexpr uint32_t binomialCoeff()
111+
{
112+
return factorial<N>() / (factorial<K>() * factorial<N - K>());
113+
}
85114
};
86115

87116
/// Helper struct for evaluating a multidimensional polynomial using compile time evaluated formula
@@ -103,7 +132,10 @@ class MultivariatePolynomialHelper : public MultivariatePolynomialParametersHelp
103132
/// evaluates the polynomial for given parameters and coordinates
104133
/// \param par parameters of the polynomials
105134
/// \param x input coordinates
106-
GPUd() static constexpr float evalPol(GPUgeneric() const float par[/*number of parameters*/], const float x[/*number of dimensions*/]) { return par[0] + loopDegrees<1>(par, x); }
135+
GPUd() static constexpr float evalPol(GPUgeneric() const float par[/*number of parameters*/], const float x[/*number of dimensions*/])
136+
{
137+
return par[0] + loopDegrees<1>(par, x);
138+
}
107139

108140
/// \return returns number of dimensions of the polynomials
109141
GPUd() static constexpr uint32_t getDim() { return Dim; }
@@ -118,19 +150,36 @@ class MultivariatePolynomialHelper : public MultivariatePolynomialParametersHelp
118150
/// computes power of 10
119151
GPUd() static constexpr uint32_t pow10(const uint32_t n) { return n == 0 ? 1 : 10 * pow10(n - 1); }
120152

153+
template <uint32_t N>
154+
GPUd() static constexpr uint32_t pow10()
155+
{
156+
if constexpr (N == 0) {
157+
return 1;
158+
} else {
159+
return 10 * pow10<N - 1>();
160+
}
161+
}
162+
121163
/// helper for modulo to extract the digit in an integer a at position b (can be obtained with pow10(digitposition)): e.g. a=1234 b=pow10(2)=100 -> returns 2
122164
GPUd() static constexpr uint32_t mod10(const uint32_t a, const uint32_t b) { return (a / b) % 10; }
123165

166+
template <uint32_t A, uint32_t B>
167+
GPUd() static constexpr uint32_t mod10()
168+
{
169+
return (A / B) % 10;
170+
}
171+
124172
/// resetting digits of pos for given position to refDigit
125173
GPUd() static constexpr uint32_t resetIndices(const uint32_t degreePol, const uint32_t pos, const uint32_t leftDigit, const uint32_t iter, const uint32_t refDigit);
126174

127-
GPUd() static constexpr uint32_t getNewPos(const uint32_t degreePol, const uint32_t pos, const uint32_t digitPos);
175+
template <uint32_t DegreePol, uint32_t Pos, uint32_t DigitPos>
176+
GPUd() static constexpr uint32_t getNewPos();
128177

129178
/// calculates term e.g. x^3*y
130179
/// \tparam DegreePol max degree of the polynomials
131180
/// \pos decoded information about the current term e.g. 1233 -> x[1]*x[2]*x[3]*x[3] (otherwise an array could be used)
132-
template <uint32_t DegreePol>
133-
GPUd() static constexpr float prodTerm(const float x[], const uint32_t pos);
181+
template <uint32_t DegreePol, uint32_t Pos>
182+
GPUd() static constexpr float prodTerm(const float x[]);
134183

135184
/// helper function for checking for interaction terms
136185
template <uint32_t DegreePol, uint32_t posNew>
@@ -203,7 +252,10 @@ class MultivariatePolynomialHelper<0, 0, false> : public MultivariatePolynomialP
203252
/// evaluating the polynomial
204253
/// \param par coefficients of the polynomial
205254
/// \param x input coordinates
206-
float evalPol(const float par[/*number of parameters*/], const float x[/*number of dimensions*/]) const { return evalPol(par, x, mDegree, mDim, mInteractionOnly); }
255+
float evalPol(const float par[/*number of parameters*/], const float x[/*number of dimensions*/]) const
256+
{
257+
return evalPol(par, x, mDegree, mDim, mInteractionOnly);
258+
}
207259

208260
/// evalutes the polynomial
209261
float evalPol(const float par[], const float x[], const uint32_t degree, const uint32_t dim, const bool interactionOnly) const;
@@ -248,35 +300,39 @@ GPUd() constexpr uint32_t MultivariatePolynomialHelper<Dim, Degree, InteractionO
248300
}
249301

250302
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
251-
GPUd() constexpr uint32_t MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::getNewPos(const uint32_t degreePol, const uint32_t pos, const uint32_t digitPos)
303+
template <uint32_t DegreePol, uint32_t Pos, uint32_t DigitPos>
304+
GPUd() constexpr uint32_t MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::getNewPos()
252305
{
253-
if (degreePol > digitPos) {
306+
if constexpr (DegreePol > DigitPos) {
254307
// check if digit of current position is at is max position
255-
if (mod10(pos, pow10(digitPos)) == Dim) {
308+
if constexpr (mod10<Pos, pow10<DigitPos>()>() == Dim) {
256309
// increase digit of left position
257-
const uint32_t leftDigit = digitPos + 1;
258-
const uint32_t posTmp = pos + pow10(leftDigit);
259-
const uint32_t refDigit = mod10(posTmp, pow10(digitPos + 1));
310+
constexpr uint32_t LeftDigit = DigitPos + 1;
311+
constexpr uint32_t PowLeftDigit = pow10<LeftDigit>();
312+
constexpr uint32_t PosTmp = Pos + PowLeftDigit;
313+
constexpr uint32_t RefDigit = mod10<PosTmp, PowLeftDigit>();
260314

261315
// resetting digits to the right if digit exceeds number of dimensions
262-
const uint32_t posReset = resetIndices(degreePol, posTmp, leftDigit - 1, degreePol - digitPos, refDigit);
316+
constexpr uint32_t PosReset = resetIndices(DegreePol, PosTmp, LeftDigit - 1, DegreePol - DigitPos, RefDigit);
263317

264318
// check next digit
265-
return getNewPos(degreePol, posReset, digitPos + 1);
319+
return getNewPos<DegreePol, PosReset, DigitPos + 1>();
320+
} else {
321+
return getNewPos<DegreePol, Pos, DigitPos + 1>();
266322
}
267-
return getNewPos(degreePol, pos, digitPos + 1);
323+
} else {
324+
return Pos;
268325
}
269-
return pos;
270326
}
271327

272328
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
273-
template <uint32_t DegreePol>
274-
GPUd() constexpr float MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::prodTerm(const float x[], const uint32_t pos)
329+
template <uint32_t DegreePol, uint32_t Pos>
330+
GPUd() constexpr float MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::prodTerm(const float x[])
275331
{
276332
if constexpr (DegreePol > 0) {
277333
// extract index of the dimension which is decoded in the digit
278-
const uint32_t index = mod10(pos, pow10(DegreePol - 1));
279-
return x[index] * prodTerm<DegreePol - 1>(x, pos);
334+
const uint32_t index = mod10<Pos, pow10<DegreePol - 1>()>();
335+
return x[index] * prodTerm<DegreePol - 1, Pos>(x);
280336
}
281337
return 1;
282338
}
@@ -286,7 +342,7 @@ template <uint32_t DegreePol, uint32_t posNew>
286342
constexpr bool MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::checkInteraction()
287343
{
288344
if constexpr (DegreePol > 1) {
289-
constexpr bool isInteraction = mod10(posNew, pow10(DegreePol - 1)) == mod10(posNew, pow10(DegreePol - 2));
345+
constexpr bool isInteraction = mod10<posNew, pow10<DegreePol - 1>()>() == mod10<posNew, pow10<DegreePol - 2>()>();
290346
if constexpr (isInteraction) {
291347
return true;
292348
}
@@ -300,16 +356,16 @@ template <uint32_t DegreePol, uint32_t Pos, uint32_t Index>
300356
GPUd() constexpr float MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::sumTerms(GPUgeneric() const float par[], const float x[])
301357
{
302358
// checking if the current position is reasonable e.g. if the max dimension is x[4]: for Pos=15 -> x[1]*x[5] the position is set to 22 -> x[2]*x[2]
303-
constexpr uint32_t posNew = getNewPos(DegreePol, Pos, 0);
304-
if constexpr (mod10(posNew, pow10(DegreePol)) != 1) {
359+
constexpr uint32_t PosNew = getNewPos<DegreePol, Pos, 0>();
360+
if constexpr (mod10<PosNew, pow10<DegreePol>()>() != 1) {
305361

306362
// check if all digits in posNew are unequal: For interaction_only terms with x[Dim]*x[Dim]... etc. can be skipped
307-
if constexpr (InteractionOnly && checkInteraction<DegreePol, posNew>()) {
308-
return sumTerms<DegreePol, posNew + 1, Index>(par, x);
363+
if constexpr (InteractionOnly && checkInteraction<DegreePol, PosNew>()) {
364+
return sumTerms<DegreePol, PosNew + 1, Index>(par, x);
365+
} else {
366+
// sum up the term for corrent term and set posotion for next combination
367+
return par[Index] * prodTerm<DegreePol, PosNew>(x) + sumTerms<DegreePol, PosNew + 1, Index + 1>(par, x);
309368
}
310-
311-
// sum up the term for corrent term and set posotion for next combination
312-
return par[Index] * prodTerm<DegreePol>(x, posNew) + sumTerms<DegreePol, posNew + 1, Index + 1>(par, x);
313369
}
314370
return 0;
315371
}
@@ -319,7 +375,7 @@ template <uint32_t DegreePol>
319375
GPUd() constexpr float MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::loopDegrees(GPUgeneric() const float par[], const float x[])
320376
{
321377
if constexpr (DegreePol <= Degree) {
322-
constexpr uint32_t index{getNParameters(DegreePol - 1, Dim, InteractionOnly)}; // offset of the index for accessing the parameters
378+
constexpr uint32_t index{getNParameters<DegreePol - 1, Dim, InteractionOnly>()}; // offset of the index for accessing the parameters
323379
return sumTerms<DegreePol, 0, index>(par, x) + loopDegrees<DegreePol + 1>(par, x);
324380
}
325381
return 0;

‎GPU/TPCFastTransformation/NDPiecewisePolynomials.h‎

Lines changed: 18 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -141,7 +141,10 @@ class NDPiecewisePolynomials : public FlatObject
141141
/// evaluate specific polynomial at given index for given coordinate
142142
/// \param x coordinates where to interpolate
143143
/// \param index index of the polynomial
144-
GPUd() float evalPol(const float x[/* Dim */], const int32_t index[/* Dim */]) const { return MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::evalPol(getParameters(index), x); }
144+
GPUd() float evalPol(const float x[/* Dim */], const int32_t index[/* Dim */]) const
145+
{
146+
return MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::evalPol(getParameters(index), x);
147+
}
145148

146149
/// \return returns min range for given dimension
147150
GPUd() float getXMin(const uint32_t dim) const { return mMin[dim]; }
@@ -215,7 +218,7 @@ class NDPiecewisePolynomials : public FlatObject
215218
#endif // !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
216219

217220
/// \return returns the total number of stored parameters
218-
uint32_t getNParameters() const { return getNPolynomials() * MultivariatePolynomialParametersHelper::getNParameters(Degree, Dim, InteractionOnly); }
221+
uint32_t getNParameters() const { return getNPolynomials() * MultivariatePolynomialParametersHelper::getNParameters<Degree, Dim, InteractionOnly>(); }
219222

220223
/// \return returns number of dimensions of the polynomials
221224
GPUd() static constexpr uint32_t getDim() { return Dim; }
@@ -241,11 +244,19 @@ class NDPiecewisePolynomials : public FlatObject
241244

242245
/// returns terms which are needed to calculate the index for the grid for given dimension
243246
/// \param dim dimension
244-
GPUd() uint32_t getTerms(const uint32_t dim) const { return (dim == 0) ? 1 : (mN[dim - 1] - 1) * getTerms(dim - 1); }
247+
template <uint32_t TermDim>
248+
GPUd() uint32_t getTerms() const
249+
{
250+
if constexpr (TermDim == 0) {
251+
return 0;
252+
} else {
253+
return (mN[TermDim - 1] - 1) * getTerms<TermDim - 1>();
254+
}
255+
}
245256

246257
/// returns index for accessing the parameter on the grid
247258
/// \param ix index per dimension
248-
GPUd() uint32_t getDataIndex(const int32_t ix[/* Dim */]) const { return getDataIndex<Dim - 1>(ix) * MultivariatePolynomialParametersHelper::getNParameters(Degree, Dim, InteractionOnly); }
259+
GPUd() uint32_t getDataIndex(const int32_t ix[/* Dim */]) const { return getDataIndex<Dim - 1>(ix) * MultivariatePolynomialParametersHelper::getNParameters<Degree, Dim, InteractionOnly>(); }
249260

250261
/// helper function to get the index
251262
template <uint32_t DimTmp>
@@ -325,7 +336,7 @@ void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::setFromContainer(cons
325336
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
326337
void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::setDefault()
327338
{
328-
const auto nParamsPerPol = MultivariatePolynomialParametersHelper::getNParameters(Degree, Dim, InteractionOnly);
339+
const auto nParamsPerPol = MultivariatePolynomialParametersHelper::getNParameters<Degree, Dim, InteractionOnly>();
329340
const auto nPols = getNPolynomials();
330341
std::vector<float> params(nParamsPerPol);
331342
params.front() = 1;
@@ -429,10 +440,10 @@ void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::setFutureBufferAddres
429440

430441
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
431442
template <uint32_t DimTmp>
432-
GPUdi() uint32_t NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::getDataIndex(const int32_t ix[/* Dim */]) const
443+
GPUd() uint32_t NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::getDataIndex(const int32_t ix[/* Dim */]) const
433444
{
434445
if constexpr (DimTmp > 0) {
435-
return ix[DimTmp] * getTerms(DimTmp) + getDataIndex<DimTmp - 1>(ix);
446+
return ix[DimTmp] * getTerms<DimTmp>() + getDataIndex<DimTmp - 1>(ix);
436447
}
437448
return ix[DimTmp];
438449
}

0 commit comments

Comments
 (0)