Line data Source code
1 : /*
2 : * File: TestOpDCT2.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on Sep 15, 2015, 4:24:44 PM
6 : *
7 : * ForBES is free software: you can redistribute it and/or modify
8 : * it under the terms of the GNU Lesser General Public License as published by
9 : * the Free Software Foundation, either version 3 of the License, or
10 : * (at your option) any later version.
11 : *
12 : * ForBES is distributed in the hope that it will be useful,
13 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : * GNU Lesser General Public License for more details.
16 : *
17 : * You should have received a copy of the GNU Lesser General Public License
18 : * along with ForBES. If not, see <http://www.gnu.org/licenses/>.
19 : */
20 :
21 : #include "TestOpDCT2.h"
22 : #include "OpAdjoint.h"
23 : #include <cmath>
24 :
25 1 : CPPUNIT_TEST_SUITE_REGISTRATION(TestOpDCT2);
26 :
27 : void testOperatorLinearity(LinearOperator* op);
28 :
29 6 : TestOpDCT2::TestOpDCT2() {
30 6 : }
31 :
32 12 : TestOpDCT2::~TestOpDCT2() {
33 12 : }
34 :
35 6 : void TestOpDCT2::setUp() {
36 6 : }
37 :
38 6 : void TestOpDCT2::tearDown() {
39 6 : }
40 :
41 1 : void TestOpDCT2::testAdj0() {
42 1 : const size_t n = 190;
43 1 : Matrix T(n, n);
44 191 : for (size_t k = 0; k < n; k++) {
45 36290 : for (size_t i = 0; i < n; i++) {
46 36100 : T.set(k, i, std::cos(static_cast<double> (k) * M_PI * (static_cast<double> (i) + 0.5) / static_cast<double> (n)));
47 : }
48 : }
49 :
50 1 : T.transpose(); // defines the adjoint of DCT-II
51 2 : Matrix z = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0);
52 :
53 2 : Matrix Tz_correct = T*z;
54 :
55 1 : LinearOperator * dct2 = new OpDCT2(n);
56 2 : Matrix Tz = dct2->callAdjoint(z);
57 1 : _ASSERT_EQ(Tz_correct, Tz);
58 :
59 2 : delete dct2;
60 1 : }
61 :
62 1 : void TestOpDCT2::testCall1() {
63 47 : for (size_t n = 14; n < 60; n++) {
64 46 : Matrix T(n, n);
65 1725 : for (size_t k = 0; k < n; k++) {
66 71070 : for (size_t i = 0; i < n; i++) {
67 69391 : T.set(k, i, std::cos(static_cast<double> (k) * M_PI * (static_cast<double> (i) + 0.5) / static_cast<double> (n)));
68 : }
69 : }
70 :
71 92 : Matrix x = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 2.0);
72 :
73 :
74 92 : Matrix y_correct = T*x;
75 :
76 46 : LinearOperator * dct2 = new OpDCT2(n);
77 92 : Matrix y = dct2->call(x);
78 :
79 92 : Matrix E = y_correct - y;
80 :
81 1725 : for (size_t j = 0; j < n; j++) {
82 1679 : _ASSERT(std::abs(E.get(j, 0)) < 1e-6);
83 : }
84 :
85 46 : delete dct2;
86 46 : }
87 1 : }
88 :
89 1 : void TestOpDCT2::testCall0() {
90 :
91 1 : const size_t n = 10;
92 :
93 1 : LinearOperator * op = new OpDCT2(n);
94 :
95 : double x_data[n] = {
96 : 0.368500000000000,
97 : -0.965700000000000,
98 : -1.332500000000000,
99 : 1.484300000000000,
100 : -0.175200000000000,
101 : -0.373400000000000,
102 : -0.512200000000000,
103 : 0.808700000000000,
104 : 0.061900000000000,
105 : 0.034600000000000
106 1 : };
107 :
108 1 : Matrix x(n, 1, x_data);
109 :
110 : double y_expected_data[n] = {
111 : -0.601000000000000,
112 : -1.162468863506907,
113 : -0.197505868217353,
114 : -0.411088627018780,
115 : 0.384982166602636,
116 : 4.028670175132235,
117 : 2.343482143524826,
118 : -0.211797156859350,
119 : -0.624017833397364,
120 : -2.578437630903463
121 1 : };
122 :
123 2 : Matrix y_expected(n, 1, y_expected_data);
124 :
125 2 : Matrix y = op->call(x);
126 2 : _ASSERT_EQ(y_expected, y);
127 :
128 :
129 1 : }
130 :
131 1 : void TestOpDCT2::testCall() {
132 1 : const size_t n = 10;
133 1 : const size_t repeat = 550;
134 1 : const double tol = 1e-8;
135 :
136 1 : LinearOperator * op = new OpDCT2(n);
137 :
138 1 : _ASSERT_EQ(n, op->dimensionIn().first);
139 1 : _ASSERT_EQ(n, op->dimensionOut().first);
140 :
141 1 : _ASSERT_NOT(op->isSelfAdjoint());
142 :
143 1 : Matrix *x = new Matrix();
144 1 : Matrix *y = new Matrix();
145 :
146 551 : for (size_t i = 0; i < repeat; i++) {
147 550 : *x = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0, Matrix::MATRIX_DENSE);
148 550 : *y = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0, Matrix::MATRIX_DENSE);
149 :
150 550 : Matrix T_y = op->call(*y);
151 1100 : Matrix Tstar_x = op->callAdjoint(*x);
152 :
153 : // make sure T_y and Tstar_x are not all zero
154 550 : bool T_y_OK = false;
155 550 : bool Tstar_x_OK = false;
156 550 : for (size_t j = 0; j < n; j++) {
157 550 : if (std::abs(T_y.get(j, 0)) > tol) {
158 550 : T_y_OK = true;
159 550 : break;
160 : }
161 : }
162 550 : for (size_t j = 0; j < n; j++) {
163 550 : if (std::abs(Tstar_x.get(j, 0)) > tol) {
164 550 : Tstar_x_OK = true;
165 550 : break;
166 : }
167 : }
168 :
169 550 : _ASSERT(T_y_OK);
170 550 : _ASSERT(Tstar_x_OK);
171 :
172 1100 : Matrix err = (*x) * T_y;
173 1100 : Matrix temp = (*y) * Tstar_x;
174 550 : Matrix::add(err, -1.0, temp, 1.0);
175 :
176 550 : double error = std::abs(err.get(0, 0));
177 550 : _ASSERT(error < tol);
178 550 : }
179 :
180 1 : delete op;
181 1 : delete x;
182 1 : delete y;
183 1 : }
184 :
185 2 : void testOperatorLinearity(LinearOperator* op) {
186 :
187 2 : double a = 10.0 * static_cast<double> (std::rand()) / static_cast<double> (RAND_MAX);
188 2 : double b = 10.0 * static_cast<double> (std::rand()) / static_cast<double> (RAND_MAX);
189 :
190 2 : Matrix x = MatrixFactory::MakeRandomMatrix(op->dimensionIn().first, op->dimensionIn().second, 0.0, 1.0);
191 4 : Matrix y = MatrixFactory::MakeRandomMatrix(op->dimensionIn().first, op->dimensionIn().second, 0.0, 1.0);
192 :
193 : // create z = ax + by
194 4 : Matrix z(x);
195 2 : Matrix::add(z, b, y, a);
196 :
197 4 : Matrix Tx = op->call(x);
198 4 : Matrix Ty = op->call(y);
199 :
200 : // create aTx + bTy
201 4 : Matrix T(Tx);
202 2 : Matrix::add(T, b, Ty, a);
203 :
204 4 : Matrix T2 = op->call(z);
205 :
206 4 : _ASSERT_EQ(T, T2);
207 :
208 2 : }
209 :
210 1 : void TestOpDCT2::testLinearity() {
211 :
212 1 : const size_t n = 50;
213 1 : LinearOperator * op = new OpDCT2(n);
214 1 : _ASSERT_EQ(n, op->dimensionIn().first);
215 1 : _ASSERT_EQ(n, op->dimensionOut().first);
216 1 : testOperatorLinearity(op);
217 1 : delete op;
218 1 : }
219 :
220 1 : void TestOpDCT2::testAdjointLinearity() {
221 1 : const size_t n = 200;
222 1 : LinearOperator * op = new OpDCT2(n);
223 1 : LinearOperator *adj = new OpAdjoint(*op);
224 1 : _ASSERT_EQ(n, adj->dimensionIn().first);
225 1 : _ASSERT_EQ(n, adj->dimensionOut().first);
226 1 : testOperatorLinearity(adj);
227 1 : delete op;
228 1 : delete adj;
229 4 : }
230 :
231 :
232 :
233 :
|
