Line data    Source code

       1             : /*
       2             :  * File:   TestOpGradient.cpp
       3             :  * Author: chung
       4             :  *
       5             :  * Created on Sep 16, 2015, 1:57:49 AM
       6             :  */
       7             : 
       8             : #include "TestOpGradient.h"
       9             : #include "OpAdjoint.h"
      10             : #include <cmath>
      11             : 
      12             : void testOperatorLinearity(LinearOperator* op);
      13             : 
      14           1 : CPPUNIT_TEST_SUITE_REGISTRATION(TestOpGradient);
      15             : 
      16           3 : TestOpGradient::TestOpGradient() {
      17           3 : }
      18             : 
      19           6 : TestOpGradient::~TestOpGradient() {
      20           6 : }
      21             : 
      22           3 : void TestOpGradient::setUp() {
      23           3 : }
      24             : 
      25           3 : void TestOpGradient::tearDown() {
      26           3 : }
      27             : 
      28           1 : void TestOpGradient::testCall() {
      29           1 :     const size_t n = 40;
      30           1 :     const double tol = 1e-8;
      31             : 
      32           1 :     LinearOperator * op = new OpGradient(n);
      33             : 
      34           1 :     _ASSERT_EQ(n, op->dimensionIn().first);
      35           1 :     _ASSERT_EQ(n - 1, op->dimensionOut().first);
      36           1 :     _ASSERT_NOT(op->isSelfAdjoint());
      37             : 
      38           1 :     Matrix x(n, 1);
      39           2 :     Matrix y(n - 1, 1);
      40             : 
      41          41 :     for (size_t i = 0; i < n; i++) {
      42          40 :         x.set(i, 0, i + 1);
      43             :     }
      44             : 
      45          40 :     for (size_t i = 0; i < n - 1; i++) {
      46          39 :         y.set(i, 0, 3 * i + 1);
      47             :     }
      48             : 
      49           2 :     Matrix Tx = op->call(x);
      50           2 :     Matrix Tstar_y = op->callAdjoint(y);
      51             : 
      52           2 :     Matrix err = y * Tx;
      53           2 :     Matrix temp = x*Tstar_y;
      54           1 :     Matrix::add(err, -1.0, temp, 1.0);
      55             : 
      56           1 :     _ASSERT(std::abs(err.get(0, 0)) < tol);
      57             : 
      58           2 :     delete op;
      59           1 : }
      60             : 
      61           2 : void testOperatorLinearity(LinearOperator* op) {
      62             : 
      63           2 :     double a = 10.0 * static_cast<double> (std::rand()) / static_cast<double> (RAND_MAX);
      64           2 :     double b = 10.0 * static_cast<double> (std::rand()) / static_cast<double> (RAND_MAX);
      65             : 
      66           2 :     Matrix x = MatrixFactory::MakeRandomMatrix(op->dimensionIn().first, op->dimensionIn().second, 0.0, 1.0);
      67           4 :     Matrix y = MatrixFactory::MakeRandomMatrix(op->dimensionIn().first, op->dimensionIn().second, 0.0, 1.0);
      68             : 
      69             :     // create z = ax + by
      70           4 :     Matrix z(x);
      71           2 :     Matrix::add(z, b, y, a);
      72             :     
      73           4 :     Matrix Tx = op->call(x);
      74           4 :     Matrix Ty = op->call(y);
      75             :     
      76             :     // create aTx + bTy
      77           4 :     Matrix T(Tx);
      78           2 :     Matrix::add(T, b, Ty, a);
      79             :     
      80           4 :     Matrix T2 = op->call(z);
      81             :     
      82           4 :     _ASSERT_EQ(T,T2);    
      83             :    
      84           2 : }
      85             : 
      86           1 : void TestOpGradient::testLinearity() {
      87           1 :     const size_t n = 50;
      88           1 :     LinearOperator * op = new OpGradient(n);
      89           1 :     testOperatorLinearity(op);
      90           1 :     delete op;
      91           1 : }
      92             : 
      93           1 : void TestOpGradient::testAdjointLinearity() {
      94           1 :     const size_t n = 50;
      95           1 :     LinearOperator * op = new OpGradient(n);
      96           1 :     LinearOperator *adj = new OpAdjoint(*op);
      97           1 :     testOperatorLinearity(adj);
      98           1 :     delete op;
      99           1 :     delete adj;
     100           4 : }