Line data Source code
1 : /**
2 : * @file test_LinearOperator.cpp
3 : *
4 : * @brief Tests for LinearOperator class
5 : *
6 : * @author Tobias Lasser - main code
7 : * @author David Frank - composite tests
8 : * @author Nikola Dinev - fixes
9 : */
10 :
11 : #include "doctest/doctest.h"
12 : #include "LinearOperator.h"
13 : #include "VolumeDescriptor.h"
14 :
15 : using namespace elsa;
16 : using namespace doctest;
17 :
18 : // mock operator, which outputs 1 for apply and 2 for applyAdjoint
19 : template <typename data_t = real_t>
20 : class MockOperator : public LinearOperator<data_t>
21 : {
22 : public:
23 : MockOperator(const DataDescriptor& domain, const DataDescriptor& range)
24 : : LinearOperator<data_t>(domain, range)
25 596 : {
26 596 : }
27 :
28 : protected:
29 : void applyImpl([[maybe_unused]] const DataContainer<data_t>& x,
30 : DataContainer<data_t>& Ax) const override
31 56 : {
32 56 : Ax = 1;
33 56 : }
34 :
35 : void applyAdjointImpl([[maybe_unused]] const DataContainer<data_t>& y,
36 : DataContainer<data_t>& Aty) const override
37 40 : {
38 40 : Aty = 3;
39 40 : }
40 :
41 : protected:
42 : MockOperator<data_t>* cloneImpl() const override
43 412 : {
44 412 : return new MockOperator<data_t>(this->getDomainDescriptor(), this->getRangeDescriptor());
45 412 : }
46 : };
47 :
48 : TYPE_TO_STRING(complex<float>);
49 : TYPE_TO_STRING(complex<double>);
50 :
51 : TEST_SUITE_BEGIN("core");
52 :
53 : TEST_CASE_TEMPLATE("LinearOperator: Testing construction", TestType, float, double, complex<float>,
54 : complex<double>)
55 12 : {
56 12 : GIVEN("DataDescriptors")
57 12 : {
58 12 : IndexVector_t numCoeff(2);
59 12 : numCoeff << 47, 11;
60 12 : IndexVector_t numCoeff2(2);
61 12 : numCoeff2 << 31, 23;
62 12 : VolumeDescriptor ddDomain(numCoeff);
63 12 : VolumeDescriptor ddRange(numCoeff2);
64 :
65 12 : WHEN("instantiating a LinearOperator")
66 12 : {
67 12 : LinearOperator<TestType> linOp(ddDomain, ddRange);
68 :
69 12 : THEN("the DataDescriptors are as expected")
70 12 : {
71 4 : REQUIRE_EQ(linOp.getDomainDescriptor(), ddDomain);
72 4 : REQUIRE_EQ(linOp.getRangeDescriptor(), ddRange);
73 4 : }
74 :
75 12 : THEN("the apply* operations throw")
76 12 : {
77 4 : DataContainer<TestType> dc(ddDomain);
78 4 : REQUIRE_THROWS_AS(linOp.apply(dc), LogicError);
79 4 : REQUIRE_THROWS_AS(linOp.applyAdjoint(dc), LogicError);
80 4 : }
81 :
82 12 : THEN("copies are good")
83 12 : {
84 4 : auto newOp = linOp;
85 4 : REQUIRE_EQ(newOp, linOp);
86 4 : }
87 12 : }
88 12 : }
89 12 : }
90 :
91 : TEST_CASE_TEMPLATE("LinearOperator: Testing clone()", TestType, float, double, complex<float>,
92 : complex<double>)
93 16 : {
94 16 : GIVEN("a LinearOperator")
95 16 : {
96 8 : IndexVector_t numCoeff(3);
97 8 : numCoeff << 23, 45, 67;
98 8 : IndexVector_t numCoeff2(2);
99 8 : numCoeff2 << 78, 90;
100 8 : VolumeDescriptor ddDomain(numCoeff);
101 8 : VolumeDescriptor ddRange(numCoeff2);
102 8 : LinearOperator<TestType> linOp(ddDomain, ddRange);
103 :
104 8 : WHEN("cloning the LinearOperator")
105 8 : {
106 8 : auto linOpClone = linOp.clone();
107 :
108 8 : THEN("everything matches")
109 8 : {
110 4 : REQUIRE_NE(linOpClone.get(), &linOp);
111 4 : REQUIRE_EQ(*linOpClone, linOp);
112 4 : }
113 :
114 8 : THEN("copies are also identical")
115 8 : {
116 4 : auto newOp = *linOpClone;
117 4 : REQUIRE_EQ(newOp, linOp);
118 4 : }
119 8 : }
120 8 : }
121 :
122 16 : GIVEN("a scalar multiplicative composite LinearOperator")
123 16 : {
124 8 : IndexVector_t numCoeff(3);
125 8 : numCoeff << 50, 41, 22;
126 8 : IndexVector_t numCoeff2(2);
127 8 : numCoeff2 << 4, 88;
128 8 : VolumeDescriptor ddDomain(numCoeff);
129 8 : VolumeDescriptor ddRange(numCoeff2);
130 8 : LinearOperator<TestType> linOp(ddDomain, ddRange);
131 8 : TestType scalar = 42;
132 :
133 8 : LinearOperator<TestType> scalarMultLinOp = scalar * linOp;
134 :
135 8 : WHEN("cloning the LinearOperator")
136 8 : {
137 8 : auto linOpClone = scalarMultLinOp.clone();
138 :
139 8 : THEN("everything matches")
140 8 : {
141 4 : REQUIRE_NE(linOpClone.get(), &scalarMultLinOp);
142 4 : REQUIRE_EQ(*linOpClone, scalarMultLinOp);
143 4 : }
144 :
145 8 : THEN("copies are also identical")
146 8 : {
147 4 : auto newOp = *linOpClone;
148 4 : REQUIRE_EQ(newOp, scalarMultLinOp);
149 4 : }
150 8 : }
151 8 : }
152 16 : }
153 :
154 : TEST_CASE_TEMPLATE("LinearOperator: Testing a leaf LinearOperator", TestType, float, double,
155 : complex<float>, complex<double>)
156 32 : {
157 32 : GIVEN("a non-adjoint leaf linear operator")
158 32 : {
159 16 : IndexVector_t numCoeff(2);
160 16 : numCoeff << 12, 23;
161 16 : IndexVector_t numCoeff2(2);
162 16 : numCoeff2 << 34, 45;
163 16 : VolumeDescriptor ddDomain(numCoeff);
164 16 : VolumeDescriptor ddRange(numCoeff2);
165 16 : MockOperator<TestType> mockOp(ddDomain, ddRange);
166 :
167 16 : auto leafOp = leaf(mockOp);
168 :
169 16 : WHEN("the operator is there")
170 16 : {
171 4 : THEN("the descriptors are set correctly")
172 4 : {
173 4 : REQUIRE_EQ(leafOp.getDomainDescriptor(), ddDomain);
174 4 : REQUIRE_EQ(leafOp.getRangeDescriptor(), ddRange);
175 4 : }
176 4 : }
177 :
178 16 : WHEN("given data")
179 16 : {
180 8 : DataContainer<TestType> dcDomain(ddDomain);
181 8 : DataContainer<TestType> dcRange(ddRange);
182 :
183 8 : THEN("the apply operations return the correct result")
184 8 : {
185 4 : auto resultApply = leafOp.apply(dcDomain);
186 6124 : for (int i = 0; i < resultApply.getSize(); ++i)
187 4 : REQUIRE_EQ(resultApply[i], static_cast<TestType>(1));
188 :
189 4 : auto resultApplyAdjoint = leafOp.applyAdjoint(dcRange);
190 1108 : for (int i = 0; i < resultApplyAdjoint.getSize(); ++i)
191 4 : REQUIRE_EQ(resultApplyAdjoint[i], static_cast<TestType>(3));
192 4 : }
193 :
194 8 : THEN("the apply operations care for appropriately sized containers")
195 8 : {
196 4 : REQUIRE_THROWS_AS(leafOp.apply(dcRange), InvalidArgumentError);
197 4 : REQUIRE_THROWS_AS(leafOp.applyAdjoint(dcDomain), InvalidArgumentError);
198 4 : }
199 8 : }
200 :
201 16 : WHEN("copying/assigning")
202 16 : {
203 4 : auto newOp = leafOp;
204 4 : auto assignedOp = leaf(newOp);
205 :
206 4 : THEN("it should be identical")
207 4 : {
208 4 : REQUIRE_EQ(newOp, leafOp);
209 4 : REQUIRE_EQ(assignedOp, leaf(newOp));
210 :
211 4 : assignedOp = newOp;
212 4 : REQUIRE_EQ(assignedOp, newOp);
213 4 : }
214 4 : }
215 16 : }
216 :
217 32 : GIVEN("an adjoint linear operator")
218 32 : {
219 16 : IndexVector_t numCoeff(2);
220 16 : numCoeff << 12, 23;
221 16 : IndexVector_t numCoeff2(2);
222 16 : numCoeff2 << 34, 45;
223 16 : VolumeDescriptor ddDomain(numCoeff);
224 16 : VolumeDescriptor ddRange(numCoeff2);
225 16 : MockOperator<TestType> mockOp(ddDomain, ddRange);
226 :
227 16 : auto adjointOp = adjoint(mockOp);
228 :
229 16 : WHEN("the operator is there")
230 16 : {
231 4 : THEN("the descriptors are set correctly")
232 4 : {
233 4 : REQUIRE_EQ(adjointOp.getDomainDescriptor(), ddRange);
234 4 : REQUIRE_EQ(adjointOp.getRangeDescriptor(), ddDomain);
235 4 : }
236 4 : }
237 :
238 16 : WHEN("given data")
239 16 : {
240 8 : DataContainer<TestType> dcDomain(ddDomain);
241 8 : DataContainer<TestType> dcRange(ddRange);
242 :
243 8 : THEN("the apply operations return the correct result")
244 8 : {
245 4 : auto resultApply = adjointOp.apply(dcRange);
246 1108 : for (int i = 0; i < resultApply.getSize(); ++i)
247 4 : REQUIRE_EQ(resultApply[i], static_cast<TestType>(3));
248 :
249 4 : auto resultApplyAdjoint = adjointOp.applyAdjoint(dcDomain);
250 6124 : for (int i = 0; i < resultApplyAdjoint.getSize(); ++i)
251 4 : REQUIRE_EQ(resultApplyAdjoint[i], static_cast<TestType>(1));
252 4 : }
253 :
254 8 : THEN("the apply operations care for appropriately sized containers")
255 8 : {
256 4 : REQUIRE_THROWS_AS(adjointOp.apply(dcDomain), InvalidArgumentError);
257 4 : REQUIRE_THROWS_AS(adjointOp.applyAdjoint(dcRange), InvalidArgumentError);
258 4 : }
259 8 : }
260 :
261 16 : WHEN("copying/assigning")
262 16 : {
263 4 : auto newOp = adjointOp;
264 4 : auto assignedOp = adjoint(newOp);
265 :
266 4 : THEN("it should be identical")
267 4 : {
268 4 : REQUIRE_EQ(newOp, adjointOp);
269 4 : REQUIRE_EQ(assignedOp, adjoint(newOp));
270 :
271 4 : assignedOp = newOp;
272 4 : REQUIRE_EQ(assignedOp, newOp);
273 4 : }
274 4 : }
275 16 : }
276 32 : }
277 :
278 : TEST_CASE_TEMPLATE("LinearOperator: Testing composite LinearOperator", TestType, float, double,
279 : complex<float>, complex<double>)
280 64 : {
281 64 : GIVEN("an additive composite linear operator")
282 64 : {
283 16 : IndexVector_t numCoeff(2);
284 16 : numCoeff << 45, 67;
285 16 : IndexVector_t numCoeff2(2);
286 16 : numCoeff2 << 13, 48;
287 16 : VolumeDescriptor ddDomain(numCoeff);
288 16 : VolumeDescriptor ddRange(numCoeff2);
289 :
290 16 : MockOperator<TestType> op1(ddDomain, ddRange);
291 16 : MockOperator<TestType> op2(ddDomain, ddRange);
292 :
293 16 : auto addOp = op1 + op2;
294 :
295 16 : WHEN("the operator is there")
296 16 : {
297 4 : THEN("the descriptors are set correctly")
298 4 : {
299 4 : REQUIRE_EQ(addOp.getDomainDescriptor(), ddDomain);
300 4 : REQUIRE_EQ(addOp.getRangeDescriptor(), ddRange);
301 4 : }
302 4 : }
303 :
304 16 : WHEN("given data")
305 16 : {
306 8 : DataContainer<TestType> dcDomain(ddDomain);
307 8 : DataContainer<TestType> dcRange(ddRange);
308 :
309 8 : THEN("the apply operations return the correct result")
310 8 : {
311 4 : auto resultApply = addOp.apply(dcDomain);
312 2500 : for (int i = 0; i < resultApply.getSize(); ++i)
313 4 : REQUIRE_EQ(resultApply[i], static_cast<TestType>(2));
314 :
315 4 : auto resultApplyAdjoint = addOp.applyAdjoint(dcRange);
316 12064 : for (int i = 0; i < resultApplyAdjoint.getSize(); ++i)
317 4 : REQUIRE_EQ(resultApplyAdjoint[i], static_cast<TestType>(6));
318 4 : }
319 :
320 8 : THEN("the apply operations care for appropriately sized containers")
321 8 : {
322 4 : REQUIRE_THROWS_AS(addOp.apply(dcRange), InvalidArgumentError);
323 4 : REQUIRE_THROWS_AS(addOp.applyAdjoint(dcDomain), InvalidArgumentError);
324 4 : }
325 8 : }
326 :
327 16 : WHEN("copying/assigning")
328 16 : {
329 4 : auto newOp = addOp;
330 4 : auto assignedOp = adjoint(newOp);
331 :
332 4 : THEN("it should be identical")
333 4 : {
334 4 : REQUIRE_EQ(newOp, addOp);
335 4 : REQUIRE_EQ(assignedOp, adjoint(newOp));
336 :
337 4 : assignedOp = newOp;
338 4 : REQUIRE_EQ(assignedOp, newOp);
339 4 : }
340 4 : }
341 16 : }
342 :
343 64 : GIVEN("a multiplicative composite linear operator")
344 64 : {
345 16 : IndexVector_t numCoeff(3);
346 16 : numCoeff << 13, 47, 69;
347 16 : IndexVector_t numCoeff2(2);
348 16 : numCoeff2 << 15, 28;
349 16 : IndexVector_t numCoeff3(4);
350 16 : numCoeff3 << 7, 30, 83, 13;
351 16 : VolumeDescriptor ddDomain(numCoeff);
352 16 : VolumeDescriptor ddMiddle(numCoeff2);
353 16 : VolumeDescriptor ddRange(numCoeff3);
354 :
355 16 : MockOperator<TestType> op1(ddDomain, ddMiddle);
356 16 : MockOperator<TestType> op2(ddMiddle, ddRange);
357 :
358 16 : auto multOp = op2 * op1;
359 :
360 16 : WHEN("the operator is there")
361 16 : {
362 4 : THEN("the descriptors are set correctly")
363 4 : {
364 4 : REQUIRE_EQ(multOp.getDomainDescriptor(), ddDomain);
365 4 : REQUIRE_EQ(multOp.getRangeDescriptor(), ddRange);
366 4 : }
367 4 : }
368 :
369 16 : WHEN("given data")
370 16 : {
371 8 : DataContainer<TestType> dcDomain(ddDomain);
372 8 : DataContainer<TestType> dcRange(ddRange);
373 :
374 8 : THEN("the apply operations return the correct result")
375 8 : {
376 4 : auto resultApply = multOp.apply(dcDomain);
377 906364 : for (int i = 0; i < resultApply.getSize(); ++i)
378 4 : REQUIRE_EQ(resultApply[i], static_cast<TestType>(1));
379 :
380 4 : auto resultApplyAdjoint = multOp.applyAdjoint(dcRange);
381 168640 : for (int i = 0; i < resultApplyAdjoint.getSize(); ++i)
382 4 : REQUIRE_EQ(resultApplyAdjoint[i], static_cast<TestType>(3));
383 4 : }
384 :
385 8 : THEN("the apply operations care for appropriately sized containers")
386 8 : {
387 4 : REQUIRE_THROWS_AS(multOp.apply(dcRange), InvalidArgumentError);
388 4 : REQUIRE_THROWS_AS(multOp.applyAdjoint(dcDomain), InvalidArgumentError);
389 4 : }
390 8 : }
391 :
392 16 : WHEN("copying/assigning")
393 16 : {
394 4 : auto newOp = multOp;
395 4 : auto assignedOp = adjoint(newOp);
396 :
397 4 : THEN("it should be identical")
398 4 : {
399 4 : REQUIRE_EQ(newOp, multOp);
400 4 : REQUIRE_EQ(assignedOp, adjoint(newOp));
401 :
402 4 : assignedOp = newOp;
403 4 : REQUIRE_EQ(assignedOp, newOp);
404 4 : }
405 4 : }
406 16 : }
407 :
408 64 : GIVEN("a scalar multiplicative composite linear operator")
409 64 : {
410 16 : IndexVector_t numCoeff(3);
411 16 : numCoeff << 13, 47, 69;
412 16 : IndexVector_t otherNumCoeff(2);
413 16 : otherNumCoeff << 15, 28;
414 16 : VolumeDescriptor ddDomain(numCoeff);
415 16 : VolumeDescriptor ddRange(otherNumCoeff);
416 :
417 16 : MockOperator op(ddDomain, ddRange);
418 16 : real_t scalar = 8;
419 :
420 16 : auto scalarMultOp = scalar * op;
421 :
422 16 : WHEN("the operator is there")
423 16 : {
424 4 : THEN("the descriptors are set correctly")
425 4 : {
426 4 : REQUIRE_EQ(scalarMultOp.getDomainDescriptor(), ddDomain);
427 4 : REQUIRE_EQ(scalarMultOp.getRangeDescriptor(), ddRange);
428 4 : }
429 4 : }
430 :
431 16 : WHEN("given data")
432 16 : {
433 8 : DataContainer dcDomain(ddDomain);
434 8 : DataContainer dcRange(ddRange);
435 :
436 8 : THEN("the apply operations return the correct result")
437 8 : {
438 4 : auto resultApply = scalarMultOp.apply(dcDomain);
439 1684 : for (int i = 0; i < resultApply.getSize(); ++i)
440 4 : REQUIRE_EQ(resultApply[i], static_cast<real_t>(8));
441 :
442 4 : auto resultApplyAdjoint = scalarMultOp.applyAdjoint(dcRange);
443 168640 : for (int i = 0; i < resultApplyAdjoint.getSize(); ++i)
444 4 : REQUIRE_EQ(resultApplyAdjoint[i], static_cast<real_t>(24));
445 4 : }
446 :
447 8 : THEN("the apply operations account for appropriately sized containers")
448 8 : {
449 4 : REQUIRE_THROWS_AS(scalarMultOp.apply(dcRange), InvalidArgumentError);
450 4 : REQUIRE_THROWS_AS(scalarMultOp.applyAdjoint(dcDomain), InvalidArgumentError);
451 4 : }
452 8 : }
453 :
454 16 : WHEN("copying/assigning")
455 16 : {
456 4 : const auto& newOp = scalarMultOp;
457 4 : auto assignedOp = adjoint(newOp);
458 :
459 4 : THEN("it should be identical")
460 4 : {
461 4 : REQUIRE_EQ(newOp, scalarMultOp);
462 4 : REQUIRE_EQ(assignedOp, adjoint(newOp));
463 :
464 4 : assignedOp = newOp;
465 4 : REQUIRE_EQ(assignedOp, newOp);
466 4 : }
467 4 : }
468 16 : }
469 :
470 64 : GIVEN("a complex composite with multiple leafs and levels")
471 64 : {
472 16 : IndexVector_t numCoeff(2);
473 16 : numCoeff << 13, 38;
474 16 : IndexVector_t numCoeff2(1);
475 16 : numCoeff2 << 16;
476 16 : IndexVector_t numCoeff3(3);
477 16 : numCoeff3 << 17, 38, 15;
478 16 : VolumeDescriptor ddDomain(numCoeff);
479 16 : VolumeDescriptor ddRange(numCoeff2);
480 16 : VolumeDescriptor ddFinalRange(numCoeff3);
481 :
482 16 : MockOperator<TestType> op1(ddDomain, ddRange);
483 16 : MockOperator<TestType> op2(ddFinalRange, ddRange);
484 16 : MockOperator<TestType> op3(ddRange, ddFinalRange);
485 :
486 16 : auto compositeOp = (op3 + adjoint(op2)) * op1;
487 :
488 16 : WHEN("the operator is there")
489 16 : {
490 4 : THEN("the descriptors are set correctly")
491 4 : {
492 4 : REQUIRE_EQ(compositeOp.getDomainDescriptor(), ddDomain);
493 4 : REQUIRE_EQ(compositeOp.getRangeDescriptor(), ddFinalRange);
494 4 : }
495 4 : }
496 :
497 16 : WHEN("given data")
498 16 : {
499 8 : DataContainer<TestType> dcDomain(ddDomain);
500 8 : DataContainer<TestType> dcFinalRange(ddFinalRange);
501 :
502 8 : THEN("the apply operations return the correct result")
503 8 : {
504 4 : auto resultApply = compositeOp.apply(dcDomain);
505 38764 : for (int i = 0; i < resultApply.getSize(); ++i)
506 4 : REQUIRE_EQ(resultApply[i], static_cast<TestType>(4));
507 :
508 4 : auto resultApplyAdjoint = compositeOp.applyAdjoint(dcFinalRange);
509 1980 : for (int i = 0; i < resultApplyAdjoint.getSize(); ++i)
510 4 : REQUIRE_EQ(resultApplyAdjoint[i], static_cast<TestType>(3));
511 4 : }
512 :
513 8 : THEN("the apply operations expect appropriately sized containers")
514 8 : {
515 4 : REQUIRE_THROWS_AS(compositeOp.apply(dcFinalRange), InvalidArgumentError);
516 4 : REQUIRE_THROWS_AS(compositeOp.applyAdjoint(dcDomain), InvalidArgumentError);
517 4 : }
518 8 : }
519 :
520 16 : WHEN("copying/assigning")
521 16 : {
522 4 : auto newOp = compositeOp;
523 4 : auto assignedOp = adjoint(newOp);
524 :
525 4 : THEN("it should be identical")
526 4 : {
527 4 : REQUIRE_EQ(newOp, compositeOp);
528 4 : REQUIRE_EQ(assignedOp, adjoint(newOp));
529 :
530 4 : assignedOp = newOp;
531 4 : REQUIRE_EQ(assignedOp, newOp);
532 4 : }
533 4 : }
534 16 : }
535 64 : }
536 :
537 : TEST_SUITE_END();
|
