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          18 : TEST_CASE_TEMPLATE("PseudoHuber: Testing without residual", TestType, float, double)
      28             : {
      29             :     using Vector = Eigen::Matrix<TestType, Eigen::Dynamic, 1>;
      30             : 
      31          12 :     GIVEN("just data (no residual)")
      32             :     {
      33          12 :         IndexVector_t numCoeff(2);
      34           6 :         numCoeff << 13, 34;
      35          12 :         VolumeDescriptor dd(numCoeff);
      36             : 
      37           6 :         real_t delta = 2;
      38             : 
      39          12 :         WHEN("instantiating")
      40             :         {
      41          12 :             PseudoHuber<TestType> func(dd, delta);
      42             : 
      43           8 :             THEN("the functional is as expected")
      44             :             {
      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             :             }
      52             : 
      53           8 :             THEN("a clone behaves as expected")
      54             :             {
      55           4 :                 auto phClone = func.clone();
      56             : 
      57           2 :                 REQUIRE_NE(phClone.get(), &func);
      58           2 :                 REQUIRE_EQ(*phClone, func);
      59             :             }
      60             : 
      61           8 :             THEN("the evaluate, gradient and Hessian work as expected")
      62             :             {
      63           4 :                 Vector dataVec(dd.getNumberOfCoefficients());
      64           2 :                 dataVec.setRandom();
      65           4 :                 DataContainer<TestType> x(dd, dataVec);
      66             : 
      67             :                 // compute the "true" values
      68           2 :                 TestType trueValue = 0;
      69           4 :                 Vector trueGrad(dd.getNumberOfCoefficients());
      70           4 :                 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             :                 }
      80             : 
      81           2 :                 REQUIRE_UNARY(checkApproxEq(func.evaluate(x), trueValue));
      82           4 :                 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             :             }
      88             :         }
      89             :     }
      90           6 : }
      91             : 
      92          18 : TEST_CASE_TEMPLATE("PseudoHuber: Testing with residual", TestType, float, double)
      93             : {
      94             :     using Vector = Eigen::Matrix<TestType, Eigen::Dynamic, 1>;
      95             : 
      96          12 :     GIVEN("a residual with data")
      97             :     {
      98             :         // linear residual
      99          12 :         IndexVector_t numCoeff(3);
     100           6 :         numCoeff << 3, 7, 11;
     101          12 :         VolumeDescriptor dd(numCoeff);
     102             : 
     103          12 :         Vector randomData(dd.getNumberOfCoefficients());
     104           6 :         randomData.setRandom();
     105          12 :         DataContainer<TestType> b(dd, randomData);
     106             : 
     107          12 :         Identity<TestType> A(dd);
     108             : 
     109          12 :         LinearResidual<TestType> linRes(A, b);
     110             : 
     111           6 :         real_t delta = static_cast<real_t>(1.5);
     112             : 
     113          12 :         WHEN("instantiating")
     114             :         {
     115          12 :             PseudoHuber<TestType> func(linRes, delta);
     116             : 
     117           8 :             THEN("the functional is as expected")
     118             :             {
     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             :             }
     125             : 
     126           8 :             THEN("a clone behaves as expected")
     127             :             {
     128           4 :                 auto phClone = func.clone();
     129             : 
     130           2 :                 REQUIRE_NE(phClone.get(), &func);
     131           2 :                 REQUIRE_EQ(*phClone, func);
     132             :             }
     133             : 
     134           8 :             THEN("the evaluate, gradient and Hessian work as expected")
     135             :             {
     136           4 :                 Vector dataVec(dd.getNumberOfCoefficients());
     137           2 :                 dataVec.setRandom();
     138           4 :                 DataContainer<TestType> x(dd, dataVec);
     139             : 
     140             :                 // compute the "true" values
     141           2 :                 TestType trueValue = 0;
     142           4 :                 Vector trueGrad(dd.getNumberOfCoefficients());
     143           4 :                 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             :                 }
     153             : 
     154           2 :                 REQUIRE_UNARY(checkApproxEq(func.evaluate(x), trueValue));
     155           4 :                 DataContainer<TestType> dcTrueGrad(dd, trueGrad);
     156           2 :                 REQUIRE_UNARY(isApprox(func.getGradient(x), dcTrueGrad));
     157             : 
     158           4 :                 auto hessian = func.getHessian(x);
     159           4 :                 auto hx = hessian.apply(x);
     160         464 :                 for (index_t i = 0; i < hx.getSize(); ++i)
     161         462 :                     REQUIRE_UNARY(checkApproxEq(hx[i], dataVec[i] * trueScale[i]));
     162             :             }
     163             :         }
     164             :     }
     165           6 : }
     166             : 
     167             : TEST_SUITE_END();