Line data    Source code

       1             : /**
       2             :  * @file test_WLSProblem.cpp
       3             :  *
       4             :  * @brief Tests for the WLSProblem class
       5             :  *
       6             :  * @author Matthias Wieczorek - initial code
       7             :  * @author Maximilian Hornung - modularization
       8             :  * @author David Frank - rewrite
       9             :  * @author Tobias Lasser - rewrite, modernization
      10             :  * @author Nikola Dinev - added tests for conversion constructor
      11             :  */
      12             : 
      13             : #include "doctest/doctest.h"
      14             : 
      15             : #include "WLSProblem.h"
      16             : #include "Identity.h"
      17             : #include "Scaling.h"
      18             : #include "L2NormPow2.h"
      19             : #include "WeightedL2NormPow2.h"
      20             : #include "Quadric.h"
      21             : #include "TikhonovProblem.h"
      22             : #include "BlockLinearOperator.h"
      23             : #include "RandomBlocksDescriptor.h"
      24             : #include "VolumeDescriptor.h"
      25             : #include "testHelpers.h"
      26             : 
      27             : using namespace elsa;
      28             : using namespace doctest;
      29             : 
      30             : TEST_SUITE_BEGIN("problems");
      31             : 
      32          26 : TEST_CASE_TEMPLATE("WLSProblems: Testing with operator and data", TestType, float, double)
      33             : {
      34          16 :     GIVEN("the operator and data")
      35             :     {
      36          16 :         IndexVector_t numCoeff(2);
      37           8 :         numCoeff << 7, 13;
      38          16 :         VolumeDescriptor dd(numCoeff);
      39             : 
      40          16 :         Eigen::Matrix<TestType, Eigen::Dynamic, 1> bVec(dd.getNumberOfCoefficients());
      41           8 :         bVec.setRandom();
      42          16 :         DataContainer<TestType> dcB(dd, bVec);
      43             : 
      44          16 :         Identity<TestType> idOp(dd);
      45             : 
      46          12 :         WHEN("setting up a ls problem without x0")
      47             :         {
      48           8 :             WLSProblem<TestType> prob(idOp, dcB);
      49             : 
      50           6 :             THEN("the clone works correctly")
      51             :             {
      52           4 :                 auto probClone = prob.clone();
      53             : 
      54           2 :                 REQUIRE_NE(probClone.get(), &prob);
      55           2 :                 REQUIRE_EQ(*probClone, prob);
      56             :             }
      57             : 
      58           6 :             THEN("the problem behaves as expected")
      59             :             {
      60           4 :                 DataContainer<TestType> dcZero(dd);
      61           2 :                 dcZero = 0;
      62           2 :                 REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), dcZero));
      63             : 
      64           2 :                 REQUIRE_UNARY(
      65             :                     checkApproxEq(prob.evaluate(), as<TestType>(0.5) * bVec.squaredNorm()));
      66           2 :                 REQUIRE_UNARY(
      67             :                     checkApproxEq(prob.getGradient(), static_cast<TestType>(-1.0f) * dcB));
      68             : 
      69           2 :                 auto hessian = prob.getHessian();
      70           2 :                 REQUIRE_UNARY(isApprox(hessian.apply(dcB), dcB));
      71             :             }
      72             :         }
      73             : 
      74          12 :         WHEN("setting up a ls problem with x0")
      75             :         {
      76           8 :             Eigen::Matrix<TestType, Eigen::Dynamic, 1> x0Vec(dd.getNumberOfCoefficients());
      77           4 :             x0Vec.setRandom();
      78           8 :             DataContainer<TestType> dcX0(dd, x0Vec);
      79             : 
      80           8 :             WLSProblem<TestType> prob(idOp, dcB, dcX0);
      81             : 
      82           6 :             THEN("the clone works correctly")
      83             :             {
      84           4 :                 auto probClone = prob.clone();
      85             : 
      86           2 :                 REQUIRE_NE(probClone.get(), &prob);
      87           2 :                 REQUIRE_EQ(*probClone, prob);
      88             :             }
      89             : 
      90           6 :             THEN("the problem behaves as expected")
      91             :             {
      92             : 
      93           2 :                 REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), dcX0));
      94             : 
      95           2 :                 REQUIRE_UNARY(checkApproxEq(prob.evaluate(),
      96             :                                             as<TestType>(0.5) * (x0Vec - bVec).squaredNorm()));
      97           2 :                 REQUIRE_UNARY(isApprox(prob.getGradient(), dcX0 - dcB));
      98             : 
      99           2 :                 auto hessian = prob.getHessian();
     100           2 :                 REQUIRE_UNARY(isApprox(hessian.apply(dcB), dcB));
     101             :             }
     102             :         }
     103             :     }
     104           8 : }
     105             : 
     106          18 : TEST_CASE_TEMPLATE("WLSProblems: Testing with weights, operator and data", TestType, float, double)
     107             : {
     108           0 :     GIVEN("weights, operator and data")
     109             :     {
     110           0 :         IndexVector_t numCoeff(3);
     111           0 :         numCoeff << 7, 13, 17;
     112           0 :         VolumeDescriptor dd(numCoeff);
     113             : 
     114           0 :         Eigen::Matrix<TestType, Eigen::Dynamic, 1> bVec(dd.getNumberOfCoefficients());
     115           0 :         bVec.setRandom();
     116           0 :         DataContainer<TestType> dcB(dd, bVec);
     117             : 
     118           0 :         Identity<TestType> idOp(dd);
     119             : 
     120           0 :         Eigen::Matrix<TestType, Eigen::Dynamic, 1> weightsVec(dd.getNumberOfCoefficients());
     121           0 :         weightsVec.setRandom();
     122           0 :         DataContainer dcWeights(dd, weightsVec);
     123           0 :         Scaling scaleOp(dd, dcWeights);
     124             : 
     125           0 :         WHEN("setting up a wls problem without x0")
     126             :         {
     127           0 :             WLSProblem<TestType> prob(scaleOp, idOp, dcB);
     128             : 
     129           0 :             THEN("the clone works correctly")
     130             :             {
     131           0 :                 auto probClone = prob.clone();
     132             : 
     133           0 :                 REQUIRE_NE(probClone.get(), &prob);
     134           0 :                 REQUIRE_EQ(*probClone, prob);
     135             :             }
     136             : 
     137           0 :             THEN("the problem behaves as expected")
     138             :             {
     139           0 :                 DataContainer<TestType> dcZero(dd);
     140           0 :                 dcZero = 0;
     141           0 :                 REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), dcZero));
     142             : 
     143           0 :                 REQUIRE_UNARY(checkApproxEq(
     144             :                     prob.evaluate(),
     145             :                     as<TestType>(0.5) * bVec.dot((weightsVec.array() * bVec.array()).matrix())));
     146           0 :                 DataContainer<TestType> tmpDc = as<TestType>(-1) * dcWeights * dcB;
     147           0 :                 REQUIRE_UNARY(isApprox(prob.getGradient(), tmpDc));
     148             : 
     149           0 :                 auto hessian = prob.getHessian();
     150           0 :                 REQUIRE_UNARY(isApprox(hessian.apply(dcB), dcWeights * dcB));
     151             :             }
     152             :         }
     153             : 
     154           0 :         WHEN("setting up a wls problem with x0")
     155             :         {
     156           0 :             Eigen::Matrix<TestType, Eigen::Dynamic, 1> x0Vec(dd.getNumberOfCoefficients());
     157           0 :             x0Vec.setRandom();
     158           0 :             DataContainer<TestType> dcX0(dd, x0Vec);
     159             : 
     160           0 :             WLSProblem<TestType> prob(scaleOp, idOp, dcB, dcX0);
     161             : 
     162           0 :             THEN("the clone works correctly")
     163             :             {
     164           0 :                 auto probClone = prob.clone();
     165             : 
     166           0 :                 REQUIRE_NE(probClone.get(), &prob);
     167           0 :                 REQUIRE_EQ(*probClone, prob);
     168             :             }
     169             : 
     170           0 :             THEN("the problem behaves as expected")
     171             :             {
     172           0 :                 DataContainer<TestType> dcZero(dd);
     173           0 :                 dcZero = 0;
     174           0 :                 REQUIRE_UNARY(isApprox(prob.getCurrentSolution(), dcX0));
     175             : 
     176           0 :                 REQUIRE_UNARY(checkApproxEq(
     177             :                     prob.evaluate(),
     178             :                     as<TestType>(0.5)
     179             :                         * (x0Vec - bVec)
     180             :                               .dot((weightsVec.array() * (x0Vec - bVec).array()).matrix())));
     181           0 :                 REQUIRE_UNARY(isApprox(prob.getGradient(), dcWeights * (dcX0 - dcB)));
     182             : 
     183           0 :                 auto hessian = prob.getHessian();
     184           0 :                 REQUIRE_UNARY(isApprox(hessian.apply(dcB), dcWeights * dcB));
     185             :             }
     186             :         }
     187             :     }
     188           0 : }
     189             : 
     190          36 : TEST_CASE_TEMPLATE("WLSProblems: Testing different optimization problems", TestType, float, double)
     191             : {
     192          20 :     GIVEN("an optimization problem with only a (w)ls data term")
     193             :     {
     194           4 :         VolumeDescriptor desc{IndexVector_t::Constant(1, 343)};
     195             : 
     196           4 :         Scaling<TestType> W{desc, static_cast<TestType>(3.0)};
     197             : 
     198           4 :         LinearResidual<TestType> residual{desc};
     199           4 :         L2NormPow2<TestType> func{residual};
     200           4 :         WeightedL2NormPow2<TestType> weightedFunc{residual, W};
     201           4 :         Problem<TestType> prob{func};
     202           4 :         Problem<TestType> weightedProb{weightedFunc};
     203             : 
     204           4 :         WHEN("converting to a wls problem")
     205             :         {
     206           4 :             WLSProblem<TestType> lsProb{prob};
     207           4 :             WLSProblem<TestType> wlsProb{weightedProb};
     208             : 
     209           4 :             THEN("only the type of the problem changes")
     210             :             {
     211           2 :                 REQUIRE_EQ(lsProb.getDataTerm(), prob.getDataTerm());
     212           2 :                 REQUIRE_EQ(lsProb.getRegularizationTerms(), prob.getRegularizationTerms());
     213           2 :                 REQUIRE_EQ(wlsProb.getDataTerm(), weightedProb.getDataTerm());
     214           2 :                 REQUIRE_EQ(wlsProb.getRegularizationTerms(), weightedProb.getRegularizationTerms());
     215             :             }
     216             :         }
     217             :     }
     218             : 
     219          20 :     GIVEN("an optimization problem with a non-(w)ls data term")
     220             :     {
     221           4 :         VolumeDescriptor desc{IndexVector_t::Constant(1, 343)};
     222             : 
     223           4 :         Quadric<TestType> quadric{desc};
     224           4 :         Problem prob{quadric};
     225             : 
     226           4 :         WHEN("converting to a WLSProblem")
     227             :         {
     228           6 :             THEN("an exception is thrown") { REQUIRE_THROWS(WLSProblem<TestType>{prob}); }
     229             :         }
     230             :     }
     231             : 
     232          20 :     GIVEN("an optimization problem with a non-(w)ls regularization term")
     233             :     {
     234           4 :         VolumeDescriptor desc{IndexVector_t::Constant(1, 343)};
     235             : 
     236           4 :         Quadric<TestType> quadric{desc};
     237           4 :         RegularizationTerm regTerm{static_cast<TestType>(5), quadric};
     238           4 :         Problem prob{L2NormPow2<TestType>{desc}, regTerm};
     239             : 
     240           4 :         WHEN("converting to a WLSProblem")
     241             :         {
     242           6 :             THEN("an exception is thrown") { REQUIRE_THROWS(WLSProblem<TestType>{prob}); }
     243             :         }
     244             :     }
     245             : 
     246          22 :     GIVEN("an optimization problem with a wls data term that has negative weighting factors")
     247             :     {
     248           8 :         VolumeDescriptor desc{IndexVector_t::Constant(1, 343)};
     249             : 
     250           8 :         Scaling<TestType> W1{desc, static_cast<TestType>(-3.0)};
     251             : 
     252           8 :         Eigen::Matrix<TestType, Eigen::Dynamic, 1> anisotropicW =
     253             :             Eigen::Matrix<TestType, Eigen::Dynamic, 1>::Constant(343, 1);
     254           4 :         anisotropicW[256] = -3.0;
     255             : 
     256           8 :         Scaling<TestType> W2{desc, DataContainer(desc, anisotropicW)};
     257             : 
     258           8 :         LinearResidual<TestType> residual{desc};
     259           8 :         WeightedL2NormPow2<TestType> weightedFunc1{residual, W1};
     260           8 :         WeightedL2NormPow2<TestType> weightedFunc2{residual, W2};
     261             : 
     262           6 :         WHEN("convering to a WLSProblem and no regularization terms are present")
     263             :         {
     264           4 :             Problem prob1{weightedFunc1};
     265           4 :             Problem prob2{weightedFunc2};
     266             : 
     267           4 :             WLSProblem converted1{prob1};
     268           4 :             WLSProblem converted2{prob2};
     269             : 
     270           4 :             THEN("only the type of the problem changes")
     271             :             {
     272           2 :                 REQUIRE_EQ(prob1.getDataTerm(), converted1.getDataTerm());
     273           2 :                 REQUIRE_EQ(prob1.getRegularizationTerms(), converted1.getRegularizationTerms());
     274           2 :                 REQUIRE_EQ(prob2.getDataTerm(), converted2.getDataTerm());
     275           2 :                 REQUIRE_EQ(prob2.getRegularizationTerms(), converted2.getRegularizationTerms());
     276             :             }
     277             :         }
     278             : 
     279           6 :         WHEN("convering to a WLSProblem and regularization terms are present")
     280             :         {
     281           4 :             RegularizationTerm regTerm{static_cast<TestType>(1.0), L2NormPow2<TestType>{desc}};
     282           4 :             Problem prob1{weightedFunc1, regTerm};
     283           4 :             Problem prob2{weightedFunc2, regTerm};
     284             : 
     285           4 :             THEN("an exception is thrown")
     286             :             {
     287           4 :                 REQUIRE_THROWS(WLSProblem{prob1});
     288           4 :                 REQUIRE_THROWS(WLSProblem{prob2});
     289             :             }
     290             :         }
     291             :     }
     292             : 
     293          20 :     GIVEN("an optimization problem with a (w)ls regularization term that has negative weighting "
     294             :           "factors")
     295             :     {
     296           4 :         VolumeDescriptor desc{IndexVector_t::Constant(1, 343)};
     297             : 
     298           4 :         Scaling<TestType> W1{desc, static_cast<TestType>(-3.0)};
     299             : 
     300           4 :         Eigen::Matrix<TestType, Eigen::Dynamic, 1> anisotropicW =
     301             :             Eigen::Matrix<TestType, Eigen::Dynamic, 1>::Constant(343, 1);
     302           2 :         anisotropicW[256] = -3.0;
     303             : 
     304           4 :         Scaling<TestType> W2{desc, DataContainer(desc, anisotropicW)};
     305             : 
     306           4 :         LinearResidual<TestType> residual{desc};
     307           4 :         WeightedL2NormPow2<TestType> weightedFunc1{residual, W1};
     308           4 :         WeightedL2NormPow2<TestType> weightedFunc2{residual, W2};
     309           4 :         L2NormPow2<TestType> nonWeightedFunc{desc};
     310             : 
     311           4 :         RegularizationTerm negWeights{static_cast<TestType>(1.0), weightedFunc1};
     312           4 :         RegularizationTerm mixedWeights{static_cast<TestType>(1.0), weightedFunc2};
     313           4 :         RegularizationTerm noWeightsNegLambda{static_cast<TestType>(-1.0), nonWeightedFunc};
     314           4 :         WHEN("convering to a WLSProblem")
     315             :         {
     316           4 :             Problem prob1{L2NormPow2<TestType>{desc}, negWeights};
     317           4 :             Problem prob2{L2NormPow2<TestType>{desc}, mixedWeights};
     318           4 :             Problem prob3{L2NormPow2<TestType>{desc}, noWeightsNegLambda};
     319             : 
     320           4 :             THEN("an exception is thrown")
     321             :             {
     322           4 :                 REQUIRE_THROWS(WLSProblem{prob1});
     323           4 :                 REQUIRE_THROWS(WLSProblem{prob2});
     324           4 :                 REQUIRE_THROWS(WLSProblem{prob3});
     325             :             }
     326             :         }
     327             :     }
     328             : 
     329          22 :     GIVEN("an OptimizationProblem with only (w)ls terms")
     330             :     {
     331           8 :         VolumeDescriptor desc{IndexVector_t::Constant(1, 343)};
     332           8 :         Eigen::Matrix<TestType, Eigen::Dynamic, 1> vec =
     333             :             Eigen::Matrix<TestType, Eigen::Dynamic, 1>::Random(343);
     334           8 :         DataContainer<TestType> b{desc, vec};
     335             : 
     336           8 :         Scaling<TestType> A{desc, static_cast<TestType>(2.0)};
     337             : 
     338           8 :         Scaling<TestType> isoW{desc, static_cast<TestType>(3.0)};
     339           8 :         Eigen::Matrix<TestType, Eigen::Dynamic, 1> vecW =
     340             :             Eigen::abs(Eigen::Matrix<TestType, Eigen::Dynamic, 1>::Random(343).array());
     341           8 :         DataContainer<TestType> dcW{desc, vecW};
     342           8 :         Scaling<TestType> nonIsoW{desc, dcW};
     343             : 
     344           8 :         std::vector<std::unique_ptr<Functional<TestType>>> dataTerms;
     345             : 
     346           4 :         dataTerms.push_back(std::make_unique<L2NormPow2<TestType>>(desc));
     347           4 :         dataTerms.push_back(std::make_unique<L2NormPow2<TestType>>(LinearResidual{b}));
     348           4 :         dataTerms.push_back(std::make_unique<L2NormPow2<TestType>>(LinearResidual{A}));
     349           4 :         dataTerms.push_back(std::make_unique<L2NormPow2<TestType>>(LinearResidual{A, b}));
     350           4 :         dataTerms.push_back(
     351             :             std::make_unique<WeightedL2NormPow2<TestType>>(LinearResidual{A, b}, isoW));
     352           4 :         dataTerms.push_back(
     353             :             std::make_unique<WeightedL2NormPow2<TestType>>(LinearResidual{A, b}, nonIsoW));
     354             : 
     355           8 :         Eigen::Matrix<TestType, Eigen::Dynamic, 1> regVec =
     356             :             Eigen::Matrix<TestType, Eigen::Dynamic, 1>::Random(343);
     357           8 :         DataContainer<TestType> bReg{desc, regVec};
     358             : 
     359           8 :         Scaling<TestType> AReg{desc, static_cast<TestType>(0.25)};
     360             : 
     361           8 :         Scaling<TestType> isoWReg{desc, static_cast<TestType>(1.5)};
     362           8 :         Eigen::Matrix<TestType, Eigen::Dynamic, 1> vecWReg =
     363             :             Eigen::abs(Eigen::Matrix<TestType, Eigen::Dynamic, 1>::Random(343).array());
     364           8 :         DataContainer<TestType> dcWReg{desc, vecWReg};
     365           8 :         Scaling<TestType> nonIsoWReg{desc, dcWReg};
     366             : 
     367           8 :         std::vector<std::unique_ptr<RegularizationTerm<TestType>>> regTerms;
     368           4 :         auto weight = static_cast<TestType>(0.5);
     369           4 :         regTerms.push_back(
     370             :             std::make_unique<RegularizationTerm<TestType>>(weight, L2NormPow2<TestType>{desc}));
     371           4 :         regTerms.push_back(std::make_unique<RegularizationTerm<TestType>>(
     372             :             weight, L2NormPow2<TestType>{LinearResidual{bReg}}));
     373           4 :         regTerms.push_back(std::make_unique<RegularizationTerm<TestType>>(
     374             :             weight, L2NormPow2<TestType>{LinearResidual{AReg}}));
     375           4 :         regTerms.push_back(std::make_unique<RegularizationTerm<TestType>>(
     376             :             weight, L2NormPow2<TestType>{LinearResidual{AReg, bReg}}));
     377           4 :         regTerms.push_back(std::make_unique<RegularizationTerm<TestType>>(
     378             :             weight, WeightedL2NormPow2{LinearResidual{AReg, bReg}, isoWReg}));
     379           4 :         regTerms.push_back(std::make_unique<RegularizationTerm<TestType>>(
     380             :             weight, WeightedL2NormPow2{LinearResidual{AReg, bReg}, nonIsoWReg}));
     381             : 
     382           4 :         std::array descriptions = {"has no operator and no vector",
     383             :                                    "has no operator, but has a vector",
     384             :                                    "has an operator, but no vector",
     385             :                                    "has an operator and a vector",
     386             :                                    "has an operator and a vector, and is weighted (isotropic)",
     387             :                                    "has an operator and a vector, and is weighted (nonisotropic)"};
     388             : 
     389          28 :         for (std::size_t i = 0; i < descriptions.size(); i++) {
     390         168 :             for (std::size_t j = 0; j < descriptions.size(); j++) {
     391         146 :                 WHEN("Trying different settings")
     392             :                 {
     393           4 :                     INFO("The data term ", descriptions[i]);
     394           4 :                     INFO("The regularization term ", descriptions[j]);
     395           4 :                     Eigen::Matrix<TestType, Eigen::Dynamic, 1> xVec =
     396             :                         Eigen::Matrix<TestType, Eigen::Dynamic, 1>::Random(343);
     397           4 :                     DataContainer<TestType> x{desc, xVec};
     398           4 :                     Problem prob{*dataTerms[i], *regTerms[j], x};
     399             : 
     400           4 :                     THEN("the problem can be converted and all operations yield the same result as "
     401             :                          "for the initial problem")
     402             :                     {
     403           4 :                         WLSProblem<TestType> converted{prob};
     404           2 :                         REQUIRE_UNARY(checkApproxEq(prob.evaluate(), converted.evaluate()));
     405             : 
     406           2 :                         auto gradDiff = prob.getGradient();
     407           2 :                         gradDiff -= converted.getGradient();
     408           2 :                         REQUIRE_UNARY(checkApproxEq(gradDiff.squaredL2Norm(), 0));
     409             :                     }
     410             :                 }
     411             :             }
     412             :         }
     413             :     }
     414             : 
     415          20 :     GIVEN("a TikhonovProblem with L2 regularization")
     416             :     {
     417           4 :         VolumeDescriptor desc{IndexVector_t::Constant(1, 343)};
     418           4 :         Eigen::Matrix<TestType, Eigen::Dynamic, 1> vec =
     419             :             Eigen::Matrix<TestType, Eigen::Dynamic, 1>::Random(343);
     420           4 :         DataContainer<TestType> b{desc, vec};
     421             : 
     422           4 :         Scaling<TestType> A{desc, static_cast<TestType>(2.0)};
     423           4 :         WLSProblem<TestType> prob{A, b};
     424             : 
     425           2 :         TestType regWeight = 4.0;
     426           6 :         TikhonovProblem<TestType> l2Reg{prob,
     427             :                                         RegularizationTerm{regWeight, L2NormPow2<TestType>{desc}}};
     428             : 
     429           4 :         THEN("the problem can be converted into a block form WLSProblem")
     430             :         {
     431           4 :             WLSProblem<TestType> conv{l2Reg};
     432             : 
     433           4 :             Scaling<TestType> lambdaScaling(desc, std::sqrt(regWeight));
     434           4 :             std::vector<std::unique_ptr<LinearOperator<TestType>>> opList(0);
     435           2 :             opList.push_back(A.clone());
     436           2 :             opList.push_back(lambdaScaling.clone());
     437           4 :             BlockLinearOperator<TestType> blockOp{opList,
     438             :                                                   BlockLinearOperator<TestType>::BlockType::ROW};
     439             : 
     440           4 :             std::vector<std::unique_ptr<DataDescriptor>> descList(0);
     441           2 :             descList.push_back(desc.clone());
     442           2 :             descList.push_back(desc.clone());
     443           4 :             RandomBlocksDescriptor vecDesc{descList};
     444           4 :             DataContainer<TestType> blockVec{vecDesc};
     445           2 :             blockVec.getBlock(0) = b;
     446           2 :             blockVec.getBlock(1) = 0;
     447             : 
     448           2 :             L2NormPow2<TestType> blockWls{LinearResidual<TestType>{blockOp, blockVec}};
     449           2 :             REQUIRE_EQ(conv.getDataTerm(), blockWls);
     450           2 :             REQUIRE_UNARY(conv.getRegularizationTerms().empty());
     451             :         }
     452             :     }
     453          18 : }
     454             : 
     455             : TEST_SUITE_END();