Line data Source code
1 : /**
2 : * @file test_CG.cpp
3 : *
4 : * @brief Tests for the CG class
5 : *
6 : * @author Nikola Dinev
7 : */
8 :
9 : #include "doctest/doctest.h"
10 :
11 : #include "CG.h"
12 : #include "Scaling.h"
13 : #include "Logger.h"
14 : #include "VolumeDescriptor.h"
15 : #include "JosephsMethod.h"
16 : #include "CircleTrajectoryGenerator.h"
17 : #include "PhantomGenerator.h"
18 : #include "TypeCasts.hpp"
19 : #include "testHelpers.h"
20 :
21 : using namespace elsa;
22 : using namespace doctest;
23 :
24 : TEST_SUITE_BEGIN("solvers");
25 :
26 : template <template <typename> typename T, typename data_t>
27 : constexpr data_t return_data_t(const T<data_t>&);
28 :
29 10 : TYPE_TO_STRING(CG<float>);
30 10 : TYPE_TO_STRING(CG<double>);
31 :
32 19 : TEST_CASE_TEMPLATE("CG: Solving a simple linear problem", TestType, CG<float>, CG<double>)
33 : {
34 : using data_t = decltype(return_data_t(std::declval<TestType>()));
35 : // eliminate the timing info from console for the tests
36 4 : Logger::setLevel(Logger::LogLevel::OFF);
37 :
38 8 : GIVEN("a linear problem")
39 : {
40 8 : IndexVector_t numCoeff(2);
41 4 : numCoeff << 13, 24;
42 8 : VolumeDescriptor dd{numCoeff};
43 :
44 8 : Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
45 4 : bVec.setRandom();
46 8 : DataContainer<data_t> dcB{dd, bVec};
47 :
48 4 : bVec.setRandom();
49 4 : bVec = bVec.cwiseAbs();
50 8 : Scaling<data_t> scalingOp{dd, DataContainer<data_t>{dd, bVec}};
51 :
52 8 : QuadricProblem<data_t> prob{scalingOp, dcB, true};
53 :
54 4 : data_t epsilon = std::numeric_limits<data_t>::epsilon();
55 :
56 6 : WHEN("setting up a cg solver")
57 : {
58 4 : TestType solver{prob, epsilon};
59 :
60 4 : THEN("the clone works correctly")
61 : {
62 4 : auto cgClone = solver.clone();
63 :
64 2 : REQUIRE_NE(cgClone.get(), &solver);
65 2 : REQUIRE_EQ(*cgClone, solver);
66 :
67 4 : AND_THEN("it works as expected")
68 : {
69 4 : auto solution = solver.solve(dd.getNumberOfCoefficients());
70 :
71 2 : DataContainer<data_t> resultsDifference = scalingOp.apply(solution) - dcB;
72 :
73 : // should have converged for the given number of iterations
74 2 : REQUIRE_LE((resultsDifference).squaredL2Norm(),
75 : epsilon * epsilon * dcB.squaredL2Norm());
76 : }
77 : }
78 : }
79 :
80 6 : WHEN("setting up a preconditioned cg solver")
81 : {
82 2 : bVec = 1 / bVec.array();
83 6 : TestType solver{prob, Scaling<data_t>{dd, DataContainer<data_t>{dd, bVec}}, epsilon};
84 :
85 4 : THEN("the clone works correctly")
86 : {
87 4 : auto cgClone = solver.clone();
88 :
89 2 : REQUIRE_NE(cgClone.get(), &solver);
90 2 : REQUIRE_EQ(*cgClone, solver);
91 :
92 4 : AND_THEN("it works as expected")
93 : {
94 : // a perfect preconditioner should allow for convergence in a single step
95 4 : auto solution = solver.solve(1);
96 :
97 2 : DataContainer<data_t> resultsDifference = scalingOp.apply(solution) - dcB;
98 :
99 : // should have converged for the given number of iterations
100 2 : REQUIRE_UNARY(checkApproxEq((resultsDifference).squaredL2Norm(),
101 : epsilon * epsilon * dcB.squaredL2Norm()));
102 : }
103 : }
104 : }
105 : }
106 4 : }
107 :
108 19 : TEST_CASE_TEMPLATE("CG: Solving a Tikhonov problem", TestType, CG<float>, CG<double>)
109 : {
110 : using data_t = decltype(return_data_t(std::declval<TestType>()));
111 : // eliminate the timing info from console for the tests
112 4 : Logger::setLevel(Logger::LogLevel::OFF);
113 :
114 8 : GIVEN("a Tikhonov problem")
115 : {
116 8 : IndexVector_t numCoeff(2);
117 4 : numCoeff << 13, 24;
118 8 : VolumeDescriptor dd{numCoeff};
119 :
120 8 : Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
121 4 : bVec.setRandom();
122 8 : DataContainer<data_t> dcB{dd, bVec};
123 :
124 4 : bVec.setRandom();
125 4 : bVec = bVec.cwiseProduct(bVec);
126 8 : Scaling<data_t> scalingOp{dd, DataContainer<data_t>{dd, bVec}};
127 :
128 4 : auto lambda = static_cast<data_t>(0.1);
129 8 : Scaling<data_t> lambdaOp{dd, lambda};
130 :
131 8 : QuadricProblem<data_t> prob{scalingOp + lambdaOp, dcB, true};
132 :
133 4 : data_t epsilon = std::numeric_limits<data_t>::epsilon();
134 :
135 6 : WHEN("setting up a cg solver")
136 : {
137 4 : TestType solver{prob, epsilon};
138 :
139 4 : THEN("the clone works correctly")
140 : {
141 4 : auto cgClone = solver.clone();
142 :
143 2 : REQUIRE_NE(cgClone.get(), &solver);
144 2 : REQUIRE_EQ(*cgClone, solver);
145 :
146 4 : AND_THEN("it works as expected")
147 : {
148 4 : auto solution = solver.solve(dd.getNumberOfCoefficients());
149 :
150 2 : DataContainer<data_t> resultsDifference =
151 : (scalingOp + lambdaOp).apply(solution) - dcB;
152 :
153 : // should have converged for the given number of iterations
154 2 : REQUIRE_LE(resultsDifference.squaredL2Norm(),
155 : epsilon * epsilon * dcB.squaredL2Norm());
156 : }
157 : }
158 : }
159 :
160 6 : WHEN("setting up a preconditioned cg solver")
161 : {
162 2 : bVec = 1 / (bVec.array() + lambda);
163 6 : TestType solver{prob, Scaling<data_t>{dd, DataContainer<data_t>{dd, bVec}}, epsilon};
164 :
165 4 : THEN("the clone works correctly")
166 : {
167 4 : auto cgClone = solver.clone();
168 :
169 2 : REQUIRE_NE(cgClone.get(), &solver);
170 2 : REQUIRE_EQ(*cgClone, solver);
171 :
172 4 : AND_THEN("it works as expected")
173 : {
174 : // a perfect preconditioner should allow for convergence in a single step
175 4 : auto solution = solver.solve(1);
176 :
177 2 : DataContainer<data_t> result = (scalingOp + lambdaOp).apply(solution);
178 :
179 : // should have converged for the given number of iterations
180 2 : REQUIRE_UNARY(checkApproxEq(result.squaredL2Norm(), dcB.squaredL2Norm()));
181 : // REQUIRE_LE(resultsDifference.squaredL2Norm(),
182 : // epsilon * epsilon * dcB.squaredL2Norm());
183 : }
184 : }
185 : }
186 : }
187 4 : }
188 :
189 1 : TEST_CASE("CG: Solving a simple phantom reconstruction")
190 : {
191 : // eliminate the timing info from console for the tests
192 1 : Logger::setLevel(Logger::LogLevel::OFF);
193 :
194 2 : GIVEN("a Phantom reconstruction problem")
195 : {
196 :
197 2 : IndexVector_t size(2);
198 1 : size << 16, 16; // TODO: determine optimal phantom size for efficient testing
199 2 : auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
200 1 : auto& volumeDescriptor = phantom.getDataDescriptor();
201 :
202 1 : index_t numAngles{30}, arc{180};
203 : auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
204 2 : numAngles, phantom.getDataDescriptor(), arc, static_cast<real_t>(size(0)) * 100.0f,
205 2 : static_cast<real_t>(size(0)));
206 :
207 2 : JosephsMethod projector(downcast<VolumeDescriptor>(volumeDescriptor), *sinoDescriptor);
208 :
209 2 : auto sinogram = projector.apply(phantom);
210 :
211 2 : WLSProblem problem(projector, sinogram);
212 1 : real_t epsilon = std::numeric_limits<real_t>::epsilon();
213 :
214 2 : WHEN("setting up a CG solver")
215 : {
216 2 : CG solver{problem, epsilon};
217 :
218 2 : THEN("the clone works correctly")
219 : {
220 2 : auto cgClone = solver.clone();
221 :
222 1 : REQUIRE_NE(cgClone.get(), &solver);
223 1 : REQUIRE_EQ(*cgClone, solver);
224 :
225 2 : AND_THEN("it works as expected")
226 : {
227 2 : auto reconstruction = solver.solve(10);
228 :
229 1 : DataContainer resultsDifference = reconstruction - phantom;
230 :
231 : // should have converged for the given number of iterations
232 : // does not converge to the optimal solution because of the regularization term
233 1 : REQUIRE_UNARY(checkApproxEq(resultsDifference.squaredL2Norm(),
234 : epsilon * epsilon * phantom.squaredL2Norm(), 0.1));
235 : }
236 : }
237 : }
238 : }
239 1 : }
240 :
241 : TEST_SUITE_END();
|