Line data    Source code

       1             : /**
       2             :  * @file test_Pseudohuber.cpp
       3             :  *
       4             :  * @brief Tests for the Huber class
       5             :  *
       6             :  * @author Matthias Wieczorek - initial code
       7             :  * @author David Frank - rewrite
       8             :  * @author Tobias Lasser - modernization
       9             :  */
      10             : 
      11             : #include <doctest/doctest.h>
      12             : 
      13             : #include "testHelpers.h"
      14             : #include "PseudoHuber.h"
      15             : #include "LinearResidual.h"
      16             : #include "LinearOperator.h"
      17             : #include "Scaling.h"
      18             : #include "Identity.h"
      19             : #include "VolumeDescriptor.h"
      20             : #include "TypeCasts.hpp"
      21             : 
      22             : using namespace elsa;
      23             : using namespace doctest;
      24             : 
      25             : TEST_SUITE_BEGIN("functionals");
      26             : 
      27             : TEST_CASE_TEMPLATE("PseudoHuber: Testing without residual", TestType, float, double)
      28           6 : {
      29           6 :     using Vector = Eigen::Matrix<TestType, Eigen::Dynamic, 1>;
      30             : 
      31           6 :     GIVEN("just data (no residual)")
      32           6 :     {
      33           6 :         IndexVector_t numCoeff(2);
      34           6 :         numCoeff << 13, 34;
      35           6 :         VolumeDescriptor dd(numCoeff);
      36             : 
      37           6 :         real_t delta = 2;
      38             : 
      39           6 :         WHEN("instantiating")
      40           6 :         {
      41           6 :             PseudoHuber<TestType> func(dd, delta);
      42             : 
      43           6 :             THEN("the functional is as expected")
      44           6 :             {
      45           2 :                 REQUIRE_EQ(func.getDomainDescriptor(), dd);
      46             : 
      47           2 :                 auto* linRes = downcast_safe<LinearResidual<TestType>>(&func.getResidual());
      48           2 :                 REQUIRE_UNARY(linRes);
      49           2 :                 REQUIRE_UNARY_FALSE(linRes->hasDataVector());
      50           2 :                 REQUIRE_UNARY_FALSE(linRes->hasOperator());
      51           2 :             }
      52             : 
      53           6 :             THEN("a clone behaves as expected")
      54           6 :             {
      55           2 :                 auto phClone = func.clone();
      56             : 
      57           2 :                 REQUIRE_NE(phClone.get(), &func);
      58           2 :                 REQUIRE_EQ(*phClone, func);
      59           2 :             }
      60             : 
      61           6 :             THEN("the evaluate, gradient and Hessian work as expected")
      62           6 :             {
      63           2 :                 Vector dataVec(dd.getNumberOfCoefficients());
      64           2 :                 dataVec.setRandom();
      65           2 :                 DataContainer<TestType> x(dd, dataVec);
      66             : 
      67             :                 // compute the "true" values
      68           2 :                 TestType trueValue = 0;
      69           2 :                 Vector trueGrad(dd.getNumberOfCoefficients());
      70           2 :                 Vector trueScale(dd.getNumberOfCoefficients());
      71         886 :                 for (index_t i = 0; i < dataVec.size(); ++i) {
      72         884 :                     TestType temp = dataVec[i] / delta;
      73         884 :                     TestType tempSq = temp * temp;
      74         884 :                     TestType sqrtOnePTempSq = std::sqrt(static_cast<TestType>(1.0) + tempSq);
      75         884 :                     trueValue += delta * delta * (sqrtOnePTempSq - static_cast<TestType>(1.0));
      76         884 :                     trueGrad[i] = dataVec[i] / sqrtOnePTempSq;
      77         884 :                     trueScale[i] = (sqrtOnePTempSq - tempSq / sqrtOnePTempSq)
      78         884 :                                    / (static_cast<TestType>(1.0) + tempSq);
      79         884 :                 }
      80             : 
      81           2 :                 REQUIRE_UNARY(checkApproxEq(func.evaluate(x), trueValue));
      82           2 :                 DataContainer<TestType> dcTrueGrad(dd, trueGrad);
      83           2 :                 REQUIRE_UNARY(isApprox(func.getGradient(x), dcTrueGrad));
      84             : 
      85           2 :                 DataContainer<TestType> dcTrueScale(dd, trueScale);
      86           2 :                 REQUIRE_EQ(func.getHessian(x), leaf(Scaling(dd, dcTrueScale)));
      87           2 :             }
      88           6 :         }
      89           6 :     }
      90           6 : }
      91             : 
      92             : TEST_CASE_TEMPLATE("PseudoHuber: Testing with residual", TestType, float, double)
      93           6 : {
      94           6 :     using Vector = Eigen::Matrix<TestType, Eigen::Dynamic, 1>;
      95             : 
      96           6 :     GIVEN("a residual with data")
      97           6 :     {
      98             :         // linear residual
      99           6 :         IndexVector_t numCoeff(3);
     100           6 :         numCoeff << 3, 7, 11;
     101           6 :         VolumeDescriptor dd(numCoeff);
     102             : 
     103           6 :         Vector randomData(dd.getNumberOfCoefficients());
     104           6 :         randomData.setRandom();
     105           6 :         DataContainer<TestType> b(dd, randomData);
     106             : 
     107           6 :         Identity<TestType> A(dd);
     108             : 
     109           6 :         LinearResidual<TestType> linRes(A, b);
     110             : 
     111           6 :         real_t delta = static_cast<real_t>(1.5);
     112             : 
     113           6 :         WHEN("instantiating")
     114           6 :         {
     115           6 :             PseudoHuber<TestType> func(linRes, delta);
     116             : 
     117           6 :             THEN("the functional is as expected")
     118           6 :             {
     119           2 :                 REQUIRE_EQ(func.getDomainDescriptor(), dd);
     120             : 
     121           2 :                 auto* lRes = downcast_safe<LinearResidual<TestType>>(&func.getResidual());
     122           2 :                 REQUIRE_UNARY(lRes);
     123           2 :                 REQUIRE_EQ(*lRes, linRes);
     124           2 :             }
     125             : 
     126           6 :             THEN("a clone behaves as expected")
     127           6 :             {
     128           2 :                 auto phClone = func.clone();
     129             : 
     130           2 :                 REQUIRE_NE(phClone.get(), &func);
     131           2 :                 REQUIRE_EQ(*phClone, func);
     132           2 :             }
     133             : 
     134           6 :             THEN("the evaluate, gradient and Hessian work as expected")
     135           6 :             {
     136           2 :                 Vector dataVec(dd.getNumberOfCoefficients());
     137           2 :                 dataVec.setRandom();
     138           2 :                 DataContainer<TestType> x(dd, dataVec);
     139             : 
     140             :                 // compute the "true" values
     141           2 :                 TestType trueValue = 0;
     142           2 :                 Vector trueGrad(dd.getNumberOfCoefficients());
     143           2 :                 Vector trueScale(dd.getNumberOfCoefficients());
     144         464 :                 for (index_t i = 0; i < dataVec.size(); ++i) {
     145         462 :                     TestType temp = (dataVec[i] - randomData[i]) / delta;
     146         462 :                     TestType tempSq = temp * temp;
     147         462 :                     TestType sqrtOnePTempSq = std::sqrt(static_cast<TestType>(1.0) + tempSq);
     148         462 :                     trueValue += delta * delta * (sqrtOnePTempSq - static_cast<TestType>(1.0));
     149         462 :                     trueGrad[i] = (dataVec[i] - randomData[i]) / sqrtOnePTempSq;
     150         462 :                     trueScale[i] = (sqrtOnePTempSq - tempSq / sqrtOnePTempSq)
     151         462 :                                    / (static_cast<TestType>(1.0) + tempSq);
     152         462 :                 }
     153             : 
     154           2 :                 REQUIRE_UNARY(checkApproxEq(func.evaluate(x), trueValue));
     155           2 :                 DataContainer<TestType> dcTrueGrad(dd, trueGrad);
     156           2 :                 REQUIRE_UNARY(isApprox(func.getGradient(x), dcTrueGrad));
     157             : 
     158           2 :                 auto hessian = func.getHessian(x);
     159           2 :                 auto hx = hessian.apply(x);
     160         464 :                 for (index_t i = 0; i < hx.getSize(); ++i)
     161           2 :                     REQUIRE_UNARY(checkApproxEq(hx[i], dataVec[i] * trueScale[i]));
     162           2 :             }
     163           6 :         }
     164           6 :     }
     165           6 : }
     166             : 
     167             : TEST_SUITE_END();