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 : TYPE_TO_STRING(CG<float>);
30 : TYPE_TO_STRING(CG<double>);
31 :
32 : TEST_CASE_TEMPLATE("CG: Solving a simple linear problem", TestType, CG<float>, CG<double>)
33 4 : {
34 4 : 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 4 : GIVEN("a linear problem")
39 4 : {
40 4 : IndexVector_t numCoeff(2);
41 4 : numCoeff << 13, 24;
42 4 : VolumeDescriptor dd{numCoeff};
43 :
44 4 : Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
45 4 : bVec.setRandom();
46 4 : DataContainer<data_t> dcB{dd, bVec};
47 :
48 4 : bVec.setRandom();
49 4 : bVec = bVec.cwiseAbs();
50 4 : Scaling<data_t> scalingOp{dd, DataContainer<data_t>{dd, bVec}};
51 :
52 4 : QuadricProblem<data_t> prob{scalingOp, dcB, true};
53 :
54 4 : data_t epsilon = std::numeric_limits<data_t>::epsilon();
55 :
56 4 : WHEN("setting up a cg solver")
57 4 : {
58 2 : TestType solver{prob, epsilon};
59 :
60 2 : THEN("the clone works correctly")
61 2 : {
62 2 : auto cgClone = solver.clone();
63 :
64 2 : REQUIRE_NE(cgClone.get(), &solver);
65 2 : REQUIRE_EQ(*cgClone, solver);
66 :
67 2 : AND_THEN("it works as expected")
68 2 : {
69 2 : 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 2 : epsilon * epsilon * dcB.squaredL2Norm());
76 2 : }
77 2 : }
78 2 : }
79 :
80 4 : WHEN("setting up a preconditioned cg solver")
81 4 : {
82 2 : bVec = 1 / bVec.array();
83 2 : TestType solver{prob, Scaling<data_t>{dd, DataContainer<data_t>{dd, bVec}}, epsilon};
84 :
85 2 : THEN("the clone works correctly")
86 2 : {
87 2 : auto cgClone = solver.clone();
88 :
89 2 : REQUIRE_NE(cgClone.get(), &solver);
90 2 : REQUIRE_EQ(*cgClone, solver);
91 :
92 2 : AND_THEN("it works as expected")
93 2 : {
94 : // a perfect preconditioner should allow for convergence in a single step
95 2 : 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 2 : epsilon * epsilon * dcB.squaredL2Norm()));
102 2 : }
103 2 : }
104 2 : }
105 4 : }
106 4 : }
107 :
108 : TEST_CASE_TEMPLATE("CG: Solving a Tikhonov problem", TestType, CG<float>, CG<double>)
109 4 : {
110 4 : 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 4 : GIVEN("a Tikhonov problem")
115 4 : {
116 4 : IndexVector_t numCoeff(2);
117 4 : numCoeff << 13, 24;
118 4 : VolumeDescriptor dd{numCoeff};
119 :
120 4 : Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
121 4 : bVec.setRandom();
122 4 : DataContainer<data_t> dcB{dd, bVec};
123 :
124 4 : bVec.setRandom();
125 4 : bVec = bVec.cwiseProduct(bVec);
126 4 : Scaling<data_t> scalingOp{dd, DataContainer<data_t>{dd, bVec}};
127 :
128 4 : auto lambda = static_cast<data_t>(0.1);
129 4 : Scaling<data_t> lambdaOp{dd, lambda};
130 :
131 4 : QuadricProblem<data_t> prob{scalingOp + lambdaOp, dcB, true};
132 :
133 4 : data_t epsilon = std::numeric_limits<data_t>::epsilon();
134 :
135 4 : WHEN("setting up a cg solver")
136 4 : {
137 2 : TestType solver{prob, epsilon};
138 :
139 2 : THEN("the clone works correctly")
140 2 : {
141 2 : auto cgClone = solver.clone();
142 :
143 2 : REQUIRE_NE(cgClone.get(), &solver);
144 2 : REQUIRE_EQ(*cgClone, solver);
145 :
146 2 : AND_THEN("it works as expected")
147 2 : {
148 2 : auto solution = solver.solve(dd.getNumberOfCoefficients());
149 :
150 2 : DataContainer<data_t> resultsDifference =
151 2 : (scalingOp + lambdaOp).apply(solution) - dcB;
152 :
153 : // should have converged for the given number of iterations
154 2 : REQUIRE_LE(resultsDifference.squaredL2Norm(),
155 2 : epsilon * epsilon * dcB.squaredL2Norm());
156 2 : }
157 2 : }
158 2 : }
159 :
160 4 : WHEN("setting up a preconditioned cg solver")
161 4 : {
162 2 : bVec = 1 / (bVec.array() + lambda);
163 2 : TestType solver{prob, Scaling<data_t>{dd, DataContainer<data_t>{dd, bVec}}, epsilon};
164 :
165 2 : THEN("the clone works correctly")
166 2 : {
167 2 : auto cgClone = solver.clone();
168 :
169 2 : REQUIRE_NE(cgClone.get(), &solver);
170 2 : REQUIRE_EQ(*cgClone, solver);
171 :
172 2 : AND_THEN("it works as expected")
173 2 : {
174 : // a perfect preconditioner should allow for convergence in a single step
175 2 : 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 2 : }
184 2 : }
185 2 : }
186 4 : }
187 4 : }
188 :
189 : TEST_CASE("CG: Solving a simple phantom reconstruction")
190 1 : {
191 : // eliminate the timing info from console for the tests
192 1 : Logger::setLevel(Logger::LogLevel::OFF);
193 :
194 1 : GIVEN("a Phantom reconstruction problem")
195 1 : {
196 :
197 1 : IndexVector_t size(2);
198 1 : size << 16, 16; // TODO: determine optimal phantom size for efficient testing
199 1 : auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
200 1 : auto& volumeDescriptor = phantom.getDataDescriptor();
201 :
202 1 : index_t numAngles{30}, arc{180};
203 1 : auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
204 1 : numAngles, phantom.getDataDescriptor(), arc, static_cast<real_t>(size(0)) * 100.0f,
205 1 : static_cast<real_t>(size(0)));
206 :
207 1 : JosephsMethod projector(downcast<VolumeDescriptor>(volumeDescriptor), *sinoDescriptor);
208 :
209 1 : auto sinogram = projector.apply(phantom);
210 :
211 1 : WLSProblem problem(projector, sinogram);
212 1 : real_t epsilon = std::numeric_limits<real_t>::epsilon();
213 :
214 1 : WHEN("setting up a CG solver")
215 1 : {
216 1 : CG solver{problem, epsilon};
217 :
218 1 : THEN("the clone works correctly")
219 1 : {
220 1 : auto cgClone = solver.clone();
221 :
222 1 : REQUIRE_NE(cgClone.get(), &solver);
223 1 : REQUIRE_EQ(*cgClone, solver);
224 :
225 1 : AND_THEN("it works as expected")
226 1 : {
227 1 : 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 1 : epsilon * epsilon * phantom.squaredL2Norm(), 0.1));
235 1 : }
236 1 : }
237 1 : }
238 1 : }
239 1 : }
240 :
241 : TEST_SUITE_END();
|
