Project
Loading...
Searching...
No Matches
testMultivarPolynomials.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
14
15#define BOOST_TEST_MODULE Test TPC Fast Transformation
16#define BOOST_TEST_MAIN
17#define BOOST_TEST_DYN_LINK
18
19#include <boost/test/unit_test.hpp>
21#include "NDPiecewisePolynomials.inc"
22#include <vector>
23
24namespace o2::gpu
25{
26
27// evaluate the polynomial of 4th degree in 5 dimensions for given coordinates and parameters
28float evalPol4_5D(const float* x, const float* par)
29{
30 return par[0] + par[1] * x[0] + par[2] * x[1] + par[3] * x[2] + par[4] * x[3] + par[5] * x[4] + par[6] * x[0] * x[0] + par[7] * x[0] * x[1] + par[8] * x[0] * x[2] + par[9] * x[0] * x[3] + par[10] * x[0] * x[4] + par[11] * x[1] * x[1] + par[12] * x[1] * x[2] + par[13] * x[1] * x[3] + par[14] * x[1] * x[4] + par[15] * x[2] * x[2] + par[16] * x[2] * x[3] + par[17] * x[2] * x[4] + par[18] * x[3] * x[3] + par[19] * x[3] * x[4] + par[20] * x[4] * x[4] + par[21] * x[0] * x[0] * x[0] + par[22] * x[0] * x[0] * x[1] + par[23] * x[0] * x[0] * x[2] + par[24] * x[0] * x[0] * x[3] + par[25] * x[0] * x[0] * x[4] + par[26] * x[0] * x[1] * x[1] + par[27] * x[0] * x[1] * x[2] + par[28] * x[0] * x[1] * x[3] + par[29] * x[0] * x[1] * x[4] + par[30] * x[0] * x[2] * x[2] + par[31] * x[0] * x[2] * x[3] + par[32] * x[0] * x[2] * x[4] + par[33] * x[0] * x[3] * x[3] + par[34] * x[0] * x[3] * x[4] + par[35] * x[0] * x[4] * x[4] + par[36] * x[1] * x[1] * x[1] + par[37] * x[1] * x[1] * x[2] + par[38] * x[1] * x[1] * x[3] + par[39] * x[1] * x[1] * x[4] + par[40] * x[1] * x[2] * x[2] + par[41] * x[1] * x[2] * x[3] + par[42] * x[1] * x[2] * x[4] + par[43] * x[1] * x[3] * x[3] + par[44] * x[1] * x[3] * x[4] + par[45] * x[1] * x[4] * x[4] + par[46] * x[2] * x[2] * x[2] + par[47] * x[2] * x[2] * x[3] + par[48] * x[2] * x[2] * x[4] + par[49] * x[2] * x[3] * x[3] + par[50] * x[2] * x[3] * x[4] + par[51] * x[2] * x[4] * x[4] + par[52] * x[3] * x[3] * x[3] + par[53] * x[3] * x[3] * x[4] + par[54] * x[3] * x[4] * x[4] + par[55] * x[4] * x[4] * x[4] + par[56] * x[0] * x[0] * x[0] * x[0] + par[57] * x[0] * x[0] * x[0] * x[1] + par[58] * x[0] * x[0] * x[0] * x[2] + par[59] * x[0] * x[0] * x[0] * x[3] + par[60] * x[0] * x[0] * x[0] * x[4] + par[61] * x[0] * x[0] * x[1] * x[1] + par[62] * x[0] * x[0] * x[1] * x[2] + par[63] * x[0] * x[0] * x[1] * x[3] + par[64] * x[0] * x[0] * x[1] * x[4] + par[65] * x[0] * x[0] * x[2] * x[2] + par[66] * x[0] * x[0] * x[2] * x[3] + par[67] * x[0] * x[0] * x[2] * x[4] + par[68] * x[0] * x[0] * x[3] * x[3] + par[69] * x[0] * x[0] * x[3] * x[4] + par[70] * x[0] * x[0] * x[4] * x[4] + par[71] * x[0] * x[1] * x[1] * x[1] + par[72] * x[0] * x[1] * x[1] * x[2] + par[73] * x[0] * x[1] * x[1] * x[3] + par[74] * x[0] * x[1] * x[1] * x[4] + par[75] * x[0] * x[1] * x[2] * x[2] + par[76] * x[0] * x[1] * x[2] * x[3] + par[77] * x[0] * x[1] * x[2] * x[4] + par[78] * x[0] * x[1] * x[3] * x[3] + par[79] * x[0] * x[1] * x[3] * x[4] + par[80] * x[0] * x[1] * x[4] * x[4] + par[81] * x[0] * x[2] * x[2] * x[2] + par[82] * x[0] * x[2] * x[2] * x[3] + par[83] * x[0] * x[2] * x[2] * x[4] + par[84] * x[0] * x[2] * x[3] * x[3] + par[85] * x[0] * x[2] * x[3] * x[4] + par[86] * x[0] * x[2] * x[4] * x[4] + par[87] * x[0] * x[3] * x[3] * x[3] + par[88] * x[0] * x[3] * x[3] * x[4] + par[89] * x[0] * x[3] * x[4] * x[4] + par[90] * x[0] * x[4] * x[4] * x[4] + par[91] * x[1] * x[1] * x[1] * x[1] + par[92] * x[1] * x[1] * x[1] * x[2] + par[93] * x[1] * x[1] * x[1] * x[3] + par[94] * x[1] * x[1] * x[1] * x[4] + par[95] * x[1] * x[1] * x[2] * x[2] + par[96] * x[1] * x[1] * x[2] * x[3] + par[97] * x[1] * x[1] * x[2] * x[4] + par[98] * x[1] * x[1] * x[3] * x[3] + par[99] * x[1] * x[1] * x[3] * x[4] + par[100] * x[1] * x[1] * x[4] * x[4] + par[101] * x[1] * x[2] * x[2] * x[2] + par[102] * x[1] * x[2] * x[2] * x[3] + par[103] * x[1] * x[2] * x[2] * x[4] + par[104] * x[1] * x[2] * x[3] * x[3] + par[105] * x[1] * x[2] * x[3] * x[4] + par[106] * x[1] * x[2] * x[4] * x[4] + par[107] * x[1] * x[3] * x[3] * x[3] + par[108] * x[1] * x[3] * x[3] * x[4] + par[109] * x[1] * x[3] * x[4] * x[4] + par[110] * x[1] * x[4] * x[4] * x[4] + par[111] * x[2] * x[2] * x[2] * x[2] + par[112] * x[2] * x[2] * x[2] * x[3] + par[113] * x[2] * x[2] * x[2] * x[4] + par[114] * x[2] * x[2] * x[3] * x[3] + par[115] * x[2] * x[2] * x[3] * x[4] + par[116] * x[2] * x[2] * x[4] * x[4] + par[117] * x[2] * x[3] * x[3] * x[3] + par[118] * x[2] * x[3] * x[3] * x[4] + par[119] * x[2] * x[3] * x[4] * x[4] + par[120] * x[2] * x[4] * x[4] * x[4] + par[121] * x[3] * x[3] * x[3] * x[3] + par[122] * x[3] * x[3] * x[3] * x[4] + par[123] * x[3] * x[3] * x[4] * x[4] + par[124] * x[3] * x[4] * x[4] * x[4] + par[125] * x[4] * x[4] * x[4] * x[4];
31}
32
33float evalPol5_5D_interactionOnly(const float* x, const float* par)
34{
35 return par[0] + x[0] * par[1] + x[1] * par[2] + x[2] * par[3] + x[3] * par[4] + x[4] * par[5] + x[0] * x[1] * par[6] + x[0] * x[2] * par[7] + x[0] * x[3] * par[8] + x[0] * x[4] * par[9] + x[1] * x[2] * par[10] + x[1] * x[3] * par[11] + x[1] * x[4] * par[12] + x[2] * x[3] * par[13] + x[2] * x[4] * par[14] + x[3] * x[4] * par[15] + x[0] * x[1] * x[2] * par[16] + x[0] * x[1] * x[3] * par[17] + x[0] * x[1] * x[4] * par[18] + x[0] * x[2] * x[3] * par[19] + x[0] * x[2] * x[4] * par[20] + x[0] * x[3] * x[4] * par[21] + x[1] * x[2] * x[3] * par[22] + x[1] * x[2] * x[4] * par[23] + x[1] * x[3] * x[4] * par[24] + x[2] * x[3] * x[4] * par[25] + x[0] * x[1] * x[2] * x[3] * par[26] + x[0] * x[1] * x[2] * x[4] * par[27] + x[0] * x[1] * x[3] * x[4] * par[28] + x[0] * x[2] * x[3] * x[4] * par[29] + x[1] * x[2] * x[3] * x[4] * par[30] + x[0] * x[1] * x[2] * x[3] * x[4] * par[31];
36}
37
38float genRand()
39{
40 const float minVal = -5;
41 const float maxVal = 5;
42 const float val = minVal + static_cast<float>(rand()) / (static_cast<float>(RAND_MAX / (maxVal - minVal)));
43 return val;
44}
45
47{
48 std::srand(std::time(nullptr));
49 const int32_t nPar5D4Deg = 126; // number of parameters
50 const int32_t nDim = 5; // dimensions
51 const int32_t nDegree = 4; // degree
52 const float abstolerance = 0.0001f; // abosulte difference between refernce to polynomial class
53
54 MultivariatePolynomial<nDim, nDegree> polCT; // compile time polynomial
55 MultivariatePolynomial<0, 0> polRT(nDim, nDegree); // run time polynomial
56
57 // compare number of parameters
58 BOOST_CHECK(nPar5D4Deg == polCT.getNParams());
59 BOOST_CHECK(nPar5D4Deg == polRT.getNParams());
60
61 float par[nPar5D4Deg]{20};
62 for (int32_t iter = 0; iter < 10; ++iter) {
63
64 // draw random parameters
65 for (int32_t i = 1; i < nPar5D4Deg; ++i) {
66 par[i] = genRand();
67 }
68
69 polCT.setParams(par);
70 polRT.setParams(par);
71
72 // compare evaluated polynomials
73 for (float a = 0; a < 1; a += 0.2f) {
74 for (float b = 0; b < 1; b += 0.2f) {
75 for (float c = 0; c < 1; c += 0.2f) {
76 for (float d = 0; d < 1; d += 0.2f) {
77 for (float e = 0; e < 1; e += 0.2f) {
78 const float arr[nDim]{a, b, c, d, e};
79 const float valCT = polCT.eval(arr);
80 const float valRT = polRT.eval(arr);
81 const float valRef = evalPol4_5D(arr, par);
82 BOOST_CHECK_SMALL(valCT - valRef, abstolerance);
83 BOOST_CHECK_SMALL(valRT - valRef, abstolerance);
84 }
85 }
86 }
87 }
88 }
89 }
90}
91
92BOOST_AUTO_TEST_CASE(Polynomials5D_InteractionOnly)
93{
94 std::srand(std::time(nullptr));
95 const int32_t nPar5D5DegInteraction = 32; // number of parameters
96 const int32_t nDim = 5; // dimensions
97 const int32_t nDegree = 5; // degree
98 const float abstolerance = 0.0001f; // abosulte difference between refernce to polynomial class
99 const bool interactionOnly = true;
100
101 MultivariatePolynomial<nDim, nDegree, interactionOnly> polCT; // compile time polynomial
102 MultivariatePolynomial<0, 0> polRT(nDim, nDegree, interactionOnly); // run time polynomial
103
104 // compare number of parameters
105 BOOST_CHECK(nPar5D5DegInteraction == polCT.getNParams());
106 BOOST_CHECK(nPar5D5DegInteraction == polRT.getNParams());
107
108 float par[nPar5D5DegInteraction]{20};
109 for (int32_t iter = 0; iter < 10; ++iter) {
110
111 // draw random parameters
112 for (int32_t i = 1; i < nPar5D5DegInteraction; ++i) {
113 par[i] = genRand();
114 }
115
116 polCT.setParams(par);
117 polRT.setParams(par);
118
119 // compare evaluated polynomials
120 for (float a = 0; a < 1; a += 0.2f) {
121 for (float b = 0; b < 1; b += 0.2f) {
122 for (float c = 0; c < 1; c += 0.2f) {
123 for (float d = 0; d < 1; d += 0.2f) {
124 for (float e = 0; e < 1; e += 0.2f) {
125 const float arr[nDim]{a, b, c, d, e};
126 const float valCT = polCT.eval(arr);
127 const float valRT = polRT.eval(arr);
128 const float valRef = evalPol5_5D_interactionOnly(arr, par);
129 BOOST_CHECK_SMALL(valCT - valRef, abstolerance);
130 BOOST_CHECK_SMALL(valRT - valRef, abstolerance);
131 }
132 }
133 }
134 }
135 }
136 }
137}
138
139BOOST_AUTO_TEST_CASE(Piecewise_polynomials)
140{
141 std::srand(std::time(nullptr));
142 const int32_t nPar5D5DegInteraction = 32; // number of parameters
143 const int32_t nDim = 5; // dimensions
144 const int32_t nDegree = 5; // degree
145 const bool interactionOnly = true; // consider only interaction terms
146
147 // reference polynomial which will be approximated by the NDPiecewisePolynomials
149
150 // define picewise polynomials with a lower degree than the original reference function
151 std::vector<float> xMin(nDim, 0);
152 std::vector<float> xMax(nDim, 1);
153 const std::vector<uint32_t> nVert(nDim, 4);
154 NDPiecewisePolynomials<nDim, nDegree - 2, interactionOnly> piecewisePol(xMin.data(), xMax.data(), nVert.data());
155
156 // define function which will be used to approximate the picewise polynomials
157 std::function<float(const double x[/* Xdim */])> refFunc = [&polCT, nDim](const double x[/* Xdim */]) {
158 float val[nDim]{};
159 std::copy(x, x + nDim, val);
160 return polCT.eval(val);
161 };
162
163 float par[nPar5D5DegInteraction]{20};
164 for (int32_t iter = 0; iter < 5; ++iter) {
165
166 // draw random parameters
167 for (int32_t i = 1; i < nPar5D5DegInteraction; ++i) {
168 par[i] = genRand();
169 }
170
171 polCT.setParams(par);
172
173 // define number of axiliary points which will be used for the fitting
174 const uint32_t nAux = 4;
175 const std::vector<uint32_t> nAxiliaryPoints(nDim, nAux);
176
177 // perform the fitting of the picewise polynomials
178 piecewisePol.performFits(refFunc, nAxiliaryPoints.data());
179
180 // compare evaluated polynomials
181 for (float a = 0; a < 1; a += 0.2f) {
182 for (float b = 0; b < 1; b += 0.2f) {
183 for (float c = 0; c < 1; c += 0.2f) {
184 for (float d = 0; d < 1; d += 0.2f) {
185 for (float e = 0; e < 1; e += 0.2f) {
186 const float arr[nDim]{a, b, c, d, e};
187 const float valCT = polCT.eval(arr);
188 const float valpiecewisePol = piecewisePol.evalUnsafe(arr);
189 if (std::abs(valCT) < 0.2) {
190 BOOST_CHECK_SMALL(valCT - valpiecewisePol, 0.02f);
191 } else {
192 BOOST_CHECK_CLOSE(valpiecewisePol, valCT, 0.5);
193 }
194 }
195 }
196 }
197 }
198 }
199
200 // check for safe evaluation outside of the grid at lower boundary
201 for (float a = xMin[0] - 10; a <= xMin[0]; a += 2.f) {
202 for (float b = xMin[1] - 10; b <= xMin[1]; b += 2.f) {
203 for (float c = xMin[2] - 10; c <= xMin[2]; c += 2.f) {
204 for (float d = xMin[3] - 10; d <= xMin[3]; d += 2.f) {
205 for (float e = xMin[4] - 10; e <= xMin[4]; e += 2.f) {
206 float arr[nDim]{a, b, c, d, e};
207 const float valUnsafeLow = piecewisePol.evalUnsafe(xMin.data());
208 const float valSafeLow = piecewisePol.eval(arr);
209 if (std::abs(valUnsafeLow) < 0.1) {
210 BOOST_CHECK_SMALL(valSafeLow - valUnsafeLow, 0.0001f);
211 } else {
212 BOOST_CHECK_CLOSE(valSafeLow, valUnsafeLow, 0.001);
213 }
214 }
215 }
216 }
217 }
218 }
219
220 // check for safe evaluation outside of the grid at upper boundary
221 const std::vector<int32_t> indexHigh(nDim, nAux - 2);
222 for (float a = xMax[0]; a <= xMax[0] + 10; a += 2.f) {
223 for (float b = xMax[1]; b <= xMax[1] + 10; b += 2.f) {
224 for (float c = xMax[2]; c <= xMax[2] + 10; c += 2.f) {
225 for (float d = xMax[3]; d <= xMax[3] + 10; d += 2.f) {
226 for (float e = xMax[4]; e <= xMax[4] + 10; e += 2.f) {
227 float arr[nDim]{a, b, c, d, e};
228 const float valUnsafeHigh = piecewisePol.evalPol(xMax.data(), indexHigh.data());
229 const float valSafeHigh = piecewisePol.eval(arr);
230 if (std::abs(valUnsafeHigh) < 0.1) {
231 BOOST_CHECK_SMALL(valSafeHigh - valUnsafeHigh, 0.0001f);
232 } else {
233 BOOST_CHECK_CLOSE(valSafeHigh, valUnsafeHigh, 0.001);
234 }
235 }
236 }
237 }
238 }
239 }
240 }
241}
242
243} // namespace o2::gpu
#define rand()
int32_t i
Timeframe settings Average distance in cm of events merged into Randomize distance around average distance of merged events Simulate correct bunch interactions instead of placing only the average number of events A value[n] sets TFlen for[n] collisions in Number of bunches per trainsort Number of bunch trains Instantaneous interaction rate Do not simulate border Total number of collisions to be placed in the interior of all time Do not select random but walk over array of events in stride steps h Event generator settings h configStandalone Use CPU for GPU Force usage of the specified GPU device no CPU fallback Continue processing after an error First event to Number of events to Number of iterations to Number of iterations to Number of initial iterations excluded from const char e
Timeframe settings Average distance in cm of events merged into Randomize distance around average distance of merged events Simulate correct bunch interactions instead of placing only the average number of events A value[n] sets TFlen for[n] collisions in Number of bunches per trainsort Number of bunch trains Instantaneous interaction rate Do not simulate border Total number of collisions to be placed in the interior of all time Do not select random but walk over array of events in stride steps h Event generator settings h configStandalone Use CPU for GPU Force usage of the specified GPU device no CPU fallback Continue processing after an error First event to Number of events to Number of iterations to Number of iterations to Number of initial iterations excluded from const char Directory with events to Events directory is and not inside events d
configStandalone Reconstruction rec_tpc multiplier for the combined cluster track error during track following area in cm for the search of for z only used if decrease removal by when above this local occupancy rowx Shift of TPC revert cluster flag edge bit distance to edge exceeds this ZOffset for secondary multiple of radius max P of electron K BB bands crossing max P of electron p BB bands crossing max P of electron d BB bands crossing max P of electron t BB bands crossing min accepted dEdxMaxTPC of the track Additive extra cluster error for Y2 if edge flag set Additive extra cluster error for Y2 if single set Additive extra cluster error for Y2 if or single set Additive extra cluster error for Z2 if or single set Extra errors to add to track when crossing depending on addErrorsCECrossing Inner portion of y range in sector that is not used in searching for extrapolated track candidates factor2YS for track merging factor2K for track merging Margin in cm of Y position when rejecting edge clusters based on uncorrected track Y Max fraction of shared clusters for tracklet Except pad from cluster finding if total number of charges in a fragment is above this Except pad from cluster finding if number of consecutive charges in a fragment is above this Threshold where a timebin is considered disabling the noisy pad check for that pad Treat signal above this minimum as saturated Smoothing factor for the exponential Highly Ionising Particle tail filter Number of timebins to use for the averaging Add additional custom track errors when crossing relative as fraction of Do N fit passes in final fit of< 0:reject at max -n hits") AddOptionRTC(rejectIFCLowRadiusCluster, uint8_t, 1, "", 0, "Reject clusters that get the IFC mask error during refit") AddOptionRTC(dEdxTruncLow, uint8_t, 2, "", 0, "Low truncation threshold, fraction of 128") AddOptionRTC(dEdxTruncHigh, uint8_t, 77, "", 0, "High truncation threshold, fraction of 128") AddOptionRTC(extrapolationTracking, int8_t, 1, "", 0, "Enable Extrapolation Tracking(prolong tracks to adjacent sectors to find short segments)") AddOptionRTC(disableRefitAttachment, uint8_t, 0, "", 0, "Bitmask to disable certain attachment steps during refit(1:attachment, 2:propagation, 4:loop following, 8:mirroring)") AddOptionRTC(rejectionStrategy, uint8_t, o2::gpu::GPUSettings::RejectionStrategyA, "", 0, "Enable rejection of TPC clusters for compression(0=no, 1=strategy A, 2=strategy B)") AddOptionRTC(mergeLoopersAfterburner, uint8_t, 1, "", 0, "Run afterburner for additional looper merging") AddOptionRTC(compressionTypeMask, uint8_t, o2::gpu::GPUSettings::CompressionFull, "", 0, "TPC Compression mode bits(1=truncate charge/width LSB, 2=differences, 4=track-model)") AddOptionRTC(compressionSortOrder, uint8_t, o2::gpu::GPUSettings::SortTime, "", 0, "Sort order of TPC compression(0=time, 1=pad, 2=Z-time-pad, 3=Z-pad-time, 4=no sorting(use incoming order))") AddOptionRTC(sigBitsCharge, uint8_t, 4, "", 0, "Number of significant bits for TPC cluster charge in compression mode 1") AddOptionRTC(sigBitsWidth, uint8_t, 3, "", 0, "Number of significant bits for TPC cluster width in compression mode 1") AddOptionRTC(dropLoopers, uint8_t, 0, "", 0, "Drop looping tracks starting from second loop") AddOptionRTC(mergerInterpolateErrors, uint8_t, 1, "", 0, "Use interpolation instead of extrapolation for chi2 based cluster rejection") AddOptionRTC(mergerInterpolateRejectAlsoOnCurrentPosition, uint8_t, 1, "", 0, "When using mergerInterpolateErrors, reject based on chi2 twice computed with interpolated and current track position starting from NDF > mergerNonInterpolateRejectMinNDF Merge tracks accross the central electrode Enable PID response Add sub threshold clusters in TPC dEdx computation OR mask of TPC flags that will reject the cluster in alternative dEdx Reject edge clusters based on uncorrected track Y during track fit Increase decrease removal by when on this number of pads from the edge h configStandalone Reconstruction rec_trd Max chi2 for TRD tracklets to be matched to a track Chi2 cut for strict matching mode Number of sigmas for ITS TPC track time error estimate Addition to search road around track prolongation along Y in cm Optimal tracklet rphi resolution in maximum chi2 per attached tracklet for TRD tracks take into account the full likelihood of the slope including exponential tails Tracks with less attached tracklets are discarded after the tracking Post trigger Pile up integration time in BCs fits all iterations but the last one in projections Use material lookup table for TPC refit Apply workarounds to avoid cyclic merge should not be needed tpc dyn configStandalone Processing proc_rtc Replace constant variables static by constexpr expressions Compile RTC in deterministic with NO_FAST_MATH flags and GPUCA_DETERMINISTIC_MODE define Use RTC to optimize GPU code h configStandalone Processing proc_rtctech Use a file lock to serialize access to the cache folder int32_t
configStandalone Reconstruction rec_tpc float
uint32_t c
Definition RawData.h:2
void setParams(const float params[])
GLint GLenum GLint x
Definition glcorearb.h:403
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLuint GLfloat * val
Definition glcorearb.h:1582
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
float evalPol5_5D_interactionOnly(const float *x, const float *par)
float evalPol4_5D(const float *x, const float *par)
GPUdoubleBinary64 b
BOOST_AUTO_TEST_CASE(Polynomials5D)
BOOST_CHECK(tree)