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