LCOV - code coverage report
Current view: top level - solvers/tests - test_CG.cpp (source / functions) Hit Total Coverage
Test: test_coverage.info.cleaned Lines: 101 101 100.0 %
Date: 2022-02-28 03:37:41 Functions: 13 13 100.0 %

          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();

Generated by: LCOV version 1.15