Line data Source code
1 : /**
2 : * @file test_ExpressionTemplates.cpp
3 : *
4 : * @brief Tests for Expression Templates
5 : *
6 : * @author Jens Petit
7 : */
8 :
9 : #define CATCH_CONFIG_ENABLE_BENCHMARKING
10 :
11 : #include "doctest/doctest.h"
12 : #include "DataContainer.h"
13 : #include "IdenticalBlocksDescriptor.h"
14 : #include <typeinfo>
15 : #include <cstdlib>
16 : #include <ctime>
17 : #include <memory>
18 : #include <cxxabi.h>
19 :
20 : using namespace elsa;
21 : using namespace doctest;
22 : static const index_t dimension = 16;
23 :
24 : // helper to print out the type
25 : template <class T>
26 0 : std::string type_name()
27 : {
28 : using typeT = typename std::remove_reference<T>::type;
29 0 : std::unique_ptr<char, void (*)(void*)> own(
30 : abi::__cxa_demangle(typeid(typeT).name(), nullptr, nullptr, nullptr), std::free);
31 :
32 0 : std::string r = own != nullptr ? own.get() : typeid(typeT).name();
33 : if (std::is_const<typeT>::value)
34 : r += " const";
35 : if (std::is_volatile<typeT>::value)
36 : r += " volatile";
37 : if (std::is_lvalue_reference<T>::value)
38 : r += "&";
39 : else if (std::is_rvalue_reference<T>::value)
40 : r += "&&";
41 0 : return r;
42 : }
43 :
44 : TEST_SUITE_BEGIN("core");
45 :
46 38 : TEST_CASE_TEMPLATE("Expression templates", TestType, float, double)
47 : {
48 44 : GIVEN("Three random data containers")
49 : {
50 12 : srand(static_cast<unsigned>(time(nullptr)));
51 :
52 24 : IndexVector_t numCoeff(3);
53 12 : numCoeff << dimension, dimension, dimension;
54 24 : VolumeDescriptor desc(numCoeff);
55 24 : DataContainer<TestType> dc(desc);
56 24 : DataContainer<TestType> dc2(desc);
57 24 : DataContainer<TestType> dc3(desc);
58 24 : DataContainer<TestType> result(desc);
59 :
60 49164 : for (index_t i = 0; i < dc.getSize(); ++i) {
61 49152 : dc[i] = static_cast<TestType>(rand()) / (static_cast<TestType>(RAND_MAX / 100.0));
62 49152 : dc2[i] = static_cast<TestType>(rand()) / (static_cast<TestType>(RAND_MAX / 100.0));
63 49152 : dc3[i] = static_cast<TestType>(rand()) / (static_cast<TestType>(RAND_MAX / 100.0));
64 : }
65 :
66 14 : WHEN("using auto with an expression")
67 : {
68 2 : auto exp = dc * dc2;
69 2 : auto exp2 = dc * dc2 - dc;
70 :
71 4 : THEN("the type is an (nested) expression type")
72 : {
73 4 : INFO(type_name<decltype(exp)>());
74 2 : INFO(type_name<decltype(exp2)>());
75 : }
76 : }
77 :
78 14 : WHEN("writing into a result DataContainer")
79 : {
80 2 : result = square(dc * dc2 - dc);
81 :
82 4 : THEN("the type is a DataContainer again") { INFO(type_name<decltype(result)>()); }
83 : }
84 :
85 14 : WHEN("Mixing expression and DataContainers")
86 : {
87 2 : Expression exp = elsa::exp(sqrt(dc) * log(dc2));
88 2 : result = dc * dc2 - dc / exp;
89 :
90 4 : THEN("the type is a DataContainer again")
91 : {
92 4 : INFO(type_name<decltype(exp)>());
93 2 : INFO(type_name<decltype(result)>());
94 : }
95 : }
96 :
97 18 : WHEN("Constructing a new DataContainer out of an expression")
98 : {
99 8 : THEN("the type is a DataContainer again")
100 : {
101 4 : DataContainer newDC = dc * dc2 + dc3 / dc2;
102 2 : INFO(type_name<decltype(newDC)>());
103 : }
104 :
105 8 : THEN("the type is a DataContainer again")
106 : {
107 4 : DataContainer newDC2 = TestType(2.8) * dc2;
108 2 : INFO(type_name<decltype(newDC2)>());
109 : }
110 :
111 8 : THEN("the type is a DataContainer again")
112 : {
113 4 : DataContainer newDC2 = dc2 * TestType(2.8);
114 2 : INFO(type_name<decltype(newDC2)>());
115 : }
116 : }
117 : }
118 :
119 44 : GIVEN("Three DataContainers")
120 : {
121 24 : IndexVector_t numCoeff(3);
122 12 : numCoeff << dimension, dimension, dimension;
123 24 : VolumeDescriptor desc(numCoeff);
124 24 : DataContainer<TestType> dc1(desc);
125 24 : DataContainer<TestType> dc2(desc);
126 24 : DataContainer<TestType> dc3(desc);
127 :
128 49164 : for (index_t i = 0; i < dc1.getSize(); ++i) {
129 49152 : dc1[i] = float(i) + 1;
130 49152 : dc2[i] = float(i) + 2;
131 49152 : dc3[i] = float(i) + 3;
132 : }
133 :
134 14 : WHEN("Performing calculations")
135 : {
136 4 : DataContainer result = dc1 + dc1 * sqrt(dc2);
137 :
138 4 : THEN("the results have to be correct")
139 : {
140 8194 : for (index_t i = 0; i < result.getSize(); ++i) {
141 8192 : REQUIRE_EQ(Approx(result[i]), dc1[i] + dc1[i] * std::sqrt(dc2[i]));
142 : }
143 : }
144 : }
145 :
146 14 : WHEN("Performing calculations")
147 : {
148 4 : DataContainer result = dc3 / dc1 * dc2 - dc3;
149 :
150 4 : THEN("the results have to be correct")
151 : {
152 8194 : for (index_t i = 0; i < result.getSize(); ++i) {
153 8192 : REQUIRE_EQ(Approx(result[i]).epsilon(0.001), dc3[i] / dc1[i] * dc2[i] - dc3[i]);
154 : }
155 : }
156 : }
157 :
158 20 : WHEN("Performing in-place calculations")
159 : {
160 10 : THEN("then the element-wise in-place addition works as expected")
161 : {
162 4 : DataContainer dc1Before = dc1;
163 2 : dc1 += dc3 * 2 / dc2;
164 8194 : for (index_t i = 0; i < dc1.getSize(); ++i) {
165 8192 : REQUIRE_EQ(Approx(dc1[i]).epsilon(0.001), dc1Before[i] + (dc3[i] * 2 / dc2[i]));
166 : }
167 : }
168 :
169 10 : THEN("then the element-wise in-place multiplication works as expected")
170 : {
171 4 : DataContainer dc1Before = dc1;
172 2 : dc1 *= dc3 * 2 / dc2;
173 8194 : for (index_t i = 0; i < dc1.getSize(); ++i) {
174 8192 : REQUIRE_EQ(Approx(dc1[i]).epsilon(0.001), dc1Before[i] * (dc3[i] * 2 / dc2[i]));
175 : }
176 : }
177 :
178 10 : THEN("then the element-wise in-place division works as expected")
179 : {
180 4 : DataContainer dc1Before = dc1;
181 2 : dc1 /= dc3 * 2 / dc2;
182 8194 : for (index_t i = 0; i < dc1.getSize(); ++i) {
183 8192 : REQUIRE_EQ(Approx(dc1[i]).epsilon(0.001), dc1Before[i] / (dc3[i] * 2 / dc2[i]));
184 : }
185 : }
186 :
187 10 : THEN("then the element-wise in-place subtraction works as expected")
188 : {
189 4 : DataContainer dc1Before = dc1;
190 2 : dc1 -= dc3 * 2 / dc2;
191 8194 : for (index_t i = 0; i < dc1.getSize(); ++i) {
192 8192 : REQUIRE_EQ(Approx(dc1[i]).epsilon(0.001), dc1Before[i] - (dc3[i] * 2 / dc2[i]));
193 : }
194 : }
195 : }
196 : }
197 :
198 38 : GIVEN("A single blocked DataContainer")
199 : {
200 12 : IndexVector_t numCoeff(2);
201 6 : numCoeff << 52, 29;
202 12 : VolumeDescriptor desc(numCoeff);
203 :
204 6 : index_t numBlocks = 7;
205 12 : IdenticalBlocksDescriptor blockDesc(numBlocks, desc);
206 :
207 12 : DataContainer<TestType> dc(blockDesc);
208 :
209 63342 : for (index_t i = 0; i < dc.getSize(); ++i) {
210 63336 : dc[i] = float(i) + 1;
211 : }
212 :
213 12 : const DataContainer constDc(blockDesc);
214 :
215 8 : WHEN("Performing regular arithmetic")
216 : {
217 4 : THEN("the calculations are correct")
218 : {
219 4 : DataContainer result = TestType(1.8) * dc + dc;
220 21114 : for (index_t i = 0; i < dc.getSize(); i++) {
221 21112 : REQUIRE_EQ(Approx(result[i]), TestType(1.8) * dc[i] + dc[i]);
222 : }
223 : }
224 : }
225 :
226 8 : WHEN("Performing blocked arithmetic")
227 : {
228 4 : THEN("the calculations are correct")
229 : {
230 4 : auto dcBlock = dc.getBlock(0);
231 4 : auto dcBlock2 = dc.getBlock(1);
232 4 : auto dcBlock3 = dc.getBlock(2);
233 :
234 4 : DataContainer result =
235 : TestType(1.8) * dcBlock + dcBlock2 / dcBlock - square(dcBlock3);
236 3018 : for (index_t i = 0; i < result.getSize(); i++) {
237 3016 : REQUIRE_EQ(Approx(result[i]), TestType(1.8) * dcBlock[i]
238 : + dcBlock2[i] / dcBlock[i]
239 : - dcBlock3[i] * dcBlock3[i]);
240 : }
241 : }
242 : }
243 :
244 8 : WHEN("Performing blocked arithmetic on a const block")
245 : {
246 4 : THEN("the calculations are correct")
247 : {
248 4 : const auto dcBlock = dc.getBlock(0);
249 :
250 4 : DataContainer result =
251 : TestType(1.8) * dcBlock + dcBlock / dcBlock - square(dcBlock);
252 3018 : for (index_t i = 0; i < result.getSize(); i++) {
253 3016 : REQUIRE_EQ(Approx(result[i]),
254 : TestType(1.8) * dc[i] + dc[i] / dc[i] - dc[i] * dc[i]);
255 : }
256 : }
257 : }
258 : }
259 :
260 34 : GIVEN("A non-blocked container")
261 : {
262 4 : IndexVector_t numCoeff(3);
263 2 : numCoeff << 52, 7, 29;
264 4 : VolumeDescriptor desc(numCoeff);
265 :
266 4 : DataContainer<TestType> dc(desc);
267 :
268 21114 : for (index_t i = 0; i < dc.getSize(); ++i) {
269 21112 : dc[i] = float(i) + 1;
270 : }
271 :
272 4 : WHEN("Creating a view and doing arithmetic")
273 : {
274 :
275 4 : IndexVector_t numCoeff(1);
276 2 : numCoeff << desc.getNumberOfCoefficients();
277 4 : VolumeDescriptor linearDesc(numCoeff);
278 4 : auto linearDc = dc.viewAs(linearDesc);
279 :
280 4 : DataContainer result =
281 : TestType(1.8) * linearDc + linearDc / linearDc - square(linearDc);
282 :
283 4 : THEN("the calculations are correct")
284 : {
285 21114 : for (index_t i = 0; i < result.getSize(); i++) {
286 21112 : REQUIRE_EQ(Approx(result[i]), TestType(1.8) * linearDc[i]
287 : + linearDc[i] / linearDc[i]
288 : - linearDc[i] * linearDc[i]);
289 : }
290 : }
291 : }
292 : }
293 :
294 : #ifdef ELSA_CUDA_VECTOR
295 : cudaDeviceReset();
296 : #endif
297 32 : }
298 :
299 : #ifdef ELSA_CUDA_VECTOR
300 : TEST_CASE_TEMPLATE("Expression templates: Testing on GPU", TestType, float, double)
301 : {
302 : GIVEN("Three random data containers")
303 : {
304 : srand(static_cast<unsigned>(time(nullptr)));
305 :
306 : IndexVector_t numCoeff(3);
307 : numCoeff << dimension, dimension, dimension;
308 : VolumeDescriptor desc(numCoeff);
309 : DataContainer<TestType> dc(desc, DataHandlerType::GPU);
310 : DataContainer<TestType> dc2(desc, DataHandlerType::GPU);
311 : DataContainer<TestType> dc3(desc, DataHandlerType::GPU);
312 : DataContainer<TestType> result(desc, DataHandlerType::GPU);
313 :
314 : for (index_t i = 0; i < dc.getSize(); ++i) {
315 : dc[i] = static_cast<TestType>(rand()) / (static_cast<TestType>(RAND_MAX / 100.0));
316 : dc2[i] = static_cast<TestType>(rand()) / (static_cast<TestType>(RAND_MAX / 100.0));
317 : dc3[i] = static_cast<TestType>(rand()) / (static_cast<TestType>(RAND_MAX / 100.0));
318 : }
319 :
320 : WHEN("using auto with an expression")
321 : {
322 : auto exp = dc * dc2;
323 : auto exp2 = dc * dc2 - dc;
324 :
325 : THEN("the type is an (nested) expression type")
326 : {
327 : INFO(type_name<decltype(exp)>());
328 : INFO(type_name<decltype(exp2)>());
329 : }
330 : }
331 :
332 : WHEN("writing into a result DataContainer")
333 : {
334 : result = square(dc * dc2 - dc);
335 :
336 : THEN("the type is a DataContainer again") { INFO(type_name<decltype(result)>()); }
337 : }
338 :
339 : WHEN("Mixing expression and DataContainers")
340 : {
341 : Expression exp = elsa::exp(sqrt(dc) * log(dc2));
342 : result = dc * dc2 - dc / exp;
343 :
344 : THEN("the type is a DataContainer again")
345 : {
346 : INFO(type_name<decltype(exp)>());
347 : INFO(type_name<decltype(result)>());
348 : }
349 : }
350 :
351 : WHEN("Constructing a new DataContainer out of an expression")
352 : {
353 : THEN("the type is a DataContainer again")
354 : {
355 : DataContainer newDC = dc * dc2 + dc3 / dc2;
356 : INFO(type_name<decltype(newDC)>());
357 : static_assert(std::is_same_v<typename decltype(newDC)::value_type, TestType>);
358 : }
359 :
360 : THEN("the type is a DataContainer again")
361 : {
362 : DataContainer newDC2 = TestType(2.8) * dc2;
363 : INFO(type_name<decltype(newDC2)>());
364 : static_assert(std::is_same_v<typename decltype(newDC2)::value_type, TestType>);
365 : }
366 :
367 : THEN("the type is a DataContainer again")
368 : {
369 : DataContainer newDC2 = dc2 * TestType(2.8);
370 : INFO(type_name<decltype(newDC2)>());
371 : static_assert(std::is_same_v<typename decltype(newDC2)::value_type, TestType>);
372 : }
373 : }
374 : }
375 :
376 : GIVEN("Three DataContainers")
377 : {
378 : IndexVector_t numCoeff(3);
379 : numCoeff << dimension, dimension, dimension;
380 : VolumeDescriptor desc(numCoeff);
381 : DataContainer dc1(desc, DataHandlerType::GPU);
382 : DataContainer dc2(desc, DataHandlerType::GPU);
383 : DataContainer dc3(desc, DataHandlerType::GPU);
384 :
385 : for (index_t i = 0; i < dc1.getSize(); ++i) {
386 : dc1[i] = float(i) + 1;
387 : dc2[i] = float(i) + 2;
388 : dc3[i] = float(i) + 3;
389 : }
390 :
391 : WHEN("Performing calculations")
392 : {
393 : DataContainer result = dc1 + dc1 * sqrt(dc2);
394 :
395 : THEN("the results have to be correct")
396 : {
397 : for (index_t i = 0; i < result.getSize(); ++i) {
398 : REQUIRE_EQ(Approx(result[i]), dc1[i] + dc1[i] * std::sqrt(dc2[i]));
399 : }
400 : }
401 : }
402 :
403 : WHEN("Performing calculations")
404 : {
405 : DataContainer result = dc3 / dc1 * dc2 - dc3;
406 :
407 : THEN("the results have to be correct")
408 : {
409 : for (index_t i = 0; i < result.getSize(); ++i) {
410 : REQUIRE_EQ(Approx(result[i]).epsilon(0.001), dc3[i] / dc1[i] * dc2[i] - dc3[i]);
411 : }
412 : }
413 : }
414 :
415 : WHEN("Performing in-place calculations")
416 : {
417 : THEN("then the element-wise in-place addition works as expected")
418 : {
419 : DataContainer dc1Before = dc1;
420 : dc1 += dc3 * 2 / dc2;
421 : for (index_t i = 0; i < dc1.getSize(); ++i) {
422 : REQUIRE_EQ(Approx(dc1[i]).epsilon(0.001), dc1Before[i] + (dc3[i] * 2 / dc2[i]));
423 : }
424 : }
425 :
426 : THEN("then the element-wise in-place multiplication works as expected")
427 : {
428 : DataContainer dc1Before = dc1;
429 : dc1 *= dc3 * 2 / dc2;
430 : for (index_t i = 0; i < dc1.getSize(); ++i) {
431 : REQUIRE_EQ(Approx(dc1[i]).epsilon(0.001), dc1Before[i] * (dc3[i] * 2 / dc2[i]));
432 : }
433 : }
434 :
435 : THEN("then the element-wise in-place division works as expected")
436 : {
437 : DataContainer dc1Before = dc1;
438 : dc1 /= dc3 * 2 / dc2;
439 : for (index_t i = 0; i < dc1.getSize(); ++i) {
440 : REQUIRE_EQ(Approx(dc1[i]).epsilon(0.001), dc1Before[i] / (dc3[i] * 2 / dc2[i]));
441 : }
442 : }
443 :
444 : THEN("then the element-wise in-place subtraction works as expected")
445 : {
446 : DataContainer dc1Before = dc1;
447 : dc1 -= dc3 * 2 / dc2;
448 : for (index_t i = 0; i < dc1.getSize(); ++i) {
449 : REQUIRE_EQ(Approx(dc1[i]).epsilon(0.001), dc1Before[i] - (dc3[i] * 2 / dc2[i]));
450 : }
451 : }
452 : }
453 : }
454 :
455 : cudaDeviceReset();
456 : }
457 : #endif
458 :
459 : TEST_SUITE_END();
|
