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