Line data Source code
1 : /*
2 : * File: TestOpComposition.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on Sep 15, 2015, 2:40:14 AM
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 "TestOpComposition.h"
22 : #include "OpReverseVector.h"
23 :
24 1 : CPPUNIT_TEST_SUITE_REGISTRATION(TestOpComposition);
25 :
26 4 : TestOpComposition::TestOpComposition() {
27 4 : }
28 :
29 8 : TestOpComposition::~TestOpComposition() {
30 8 : }
31 :
32 4 : void TestOpComposition::setUp() {
33 4 : }
34 :
35 4 : void TestOpComposition::tearDown() {
36 4 : }
37 :
38 1 : void TestOpComposition::testCall() {
39 1 : const size_t n = 30;
40 1 : const double tol = 1e-8;
41 1 : Matrix A = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_DENSE);
42 2 : Matrix B = MatrixFactory::MakeRandomMatrix(n, n, 0.0, 1.0, Matrix::MATRIX_DENSE);
43 1 : MatrixOperator * T1 = new MatrixOperator(A); // T1(x) = A*x
44 1 : MatrixOperator * T2 = new MatrixOperator(B); // T2(x) = B*x
45 :
46 1 : OpComposition *G = new OpComposition(*T1, *T2); // G(x) = T1(T2(x))
47 :
48 2 : Matrix x = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 1.0, Matrix::MATRIX_DENSE); // random vector x
49 :
50 2 : Matrix z = G -> call(x); // z = G(x) = T1(T2(x))
51 2 : Matrix t = T2 -> call(x); // t = T2(x)
52 2 : Matrix u = T1 -> call(t);
53 2 : Matrix err = z - u;
54 31 : for (size_t i = 0; i < n; i++) {
55 30 : _ASSERT_NUM_EQ(err.get(i, 0), 0.0, tol);
56 : }
57 :
58 : /* --- and now with different sizes */
59 : // x: m-by-1
60 : // A: p-by-n
61 : // B: n-ny-m
62 1 : const size_t m = 8;
63 1 : const size_t p = 15;
64 1 : A = MatrixFactory::MakeRandomMatrix(p, n, 0.0, 1.0, Matrix::MATRIX_DENSE);
65 1 : B = MatrixFactory::MakeRandomMatrix(n, m, 0.0, 1.0, Matrix::MATRIX_DENSE);
66 :
67 1 : T1 = new MatrixOperator(A);
68 1 : T2 = new MatrixOperator(B);
69 :
70 :
71 1 : G = new OpComposition(*T1, *T2);
72 :
73 1 : _ASSERT_EQ(p, T1->dimensionOut().first);
74 1 : _ASSERT_EQ(static_cast<size_t> (1), T1->dimensionOut().second);
75 :
76 1 : _ASSERT_EQ(p, G->dimensionOut().first);
77 1 : _ASSERT_EQ(static_cast<size_t> (1), G->dimensionOut().second);
78 :
79 1 : _ASSERT_EQ(m, G->dimensionIn().first);
80 1 : _ASSERT_EQ(static_cast<size_t> (1), G->dimensionIn().second);
81 :
82 1 : x = MatrixFactory::MakeRandomMatrix(m, 1, 0.0, 1.0, Matrix::MATRIX_DENSE);
83 :
84 1 : z = G->call(x);
85 1 : t = T2 -> call(x); // t = T2(x)
86 1 : u = T1 -> call(t);
87 1 : err = z - u;
88 16 : for (size_t i = 0; i < G->dimensionOut().first; i++) {
89 15 : _ASSERT_NUM_EQ(err.get(i, 0), 0.0, tol);
90 : }
91 :
92 1 : delete T1;
93 1 : delete T2;
94 2 : delete G;
95 :
96 1 : }
97 :
98 1 : void TestOpComposition::testCall2() {
99 1 : size_t n = 10;
100 1 : size_t m = 15;
101 :
102 1 : LinearOperator * rev_op = new OpReverseVector(n);
103 1 : Matrix A = MatrixFactory::MakeRandomMatrix(m, n, -1.0, 2.0, Matrix::MATRIX_DENSE);
104 1 : LinearOperator * mat_op = new MatrixOperator(A);
105 :
106 1 : LinearOperator * op = new OpComposition(*mat_op, *rev_op);
107 :
108 2 : Matrix y = MatrixFactory::MakeRandomMatrix(m, 1, 0.0, 3.0, Matrix::MATRIX_DENSE);
109 :
110 2 : Matrix x = MatrixFactory::MakeRandomMatrix(n, 1, 0.0, 3.0, Matrix::MATRIX_DENSE);
111 :
112 1 : double alpha = -3.52540;
113 1 : double gamma = 1.41451;
114 2 : Matrix y_copy(y);
115 1 : int status = op->call(y, alpha, x, gamma);
116 1 : _ASSERT_EQ(ForBESUtils::STATUS_OK, status);
117 :
118 : // verify result;
119 2 : Matrix x_rev = rev_op->call(x); // rev(x)
120 2 : Matrix t = A*x_rev; //
121 1 : t *= alpha;
122 1 : y_copy *= gamma;
123 1 : y_copy += t;
124 1 : _ASSERT_EQ(y_copy, y);
125 :
126 1 : delete rev_op;
127 1 : delete mat_op;
128 2 : delete op;
129 1 : }
130 :
131 1 : void TestOpComposition::testCallAdjoint() {
132 1 : const size_t p = 8;
133 1 : const size_t n = 15;
134 1 : const size_t m = 6;
135 :
136 1 : Matrix A = MatrixFactory::MakeRandomMatrix(p, n, 0.0, 1.0, Matrix::MATRIX_DENSE);
137 2 : Matrix B = MatrixFactory::MakeRandomMatrix(n, m, 0.0, 1.0, Matrix::MATRIX_DENSE);
138 :
139 :
140 1 : LinearOperator * A_op = new MatrixOperator(A);
141 1 : LinearOperator * B_op = new MatrixOperator(B);
142 :
143 1 : LinearOperator * AB_op = new OpComposition(*A_op, *B_op);
144 :
145 2 : Matrix x = MatrixFactory::MakeRandomMatrix(p, 1, 0.0, 1.0, Matrix::MATRIX_DENSE);
146 :
147 2 : Matrix AB_op_adj_x = AB_op->callAdjoint(x);
148 :
149 2 : Matrix AB_op_adj_x_correct;
150 1 : B.transpose();
151 1 : A.transpose();
152 1 : AB_op_adj_x_correct = B * A * x;
153 1 : _ASSERT_EQ(AB_op_adj_x_correct, AB_op_adj_x);
154 1 : _ASSERT_NOT(AB_op->isSelfAdjoint());
155 1 : delete A_op;
156 1 : delete B_op;
157 2 : delete AB_op;
158 :
159 1 : }
160 :
161 1 : void TestOpComposition::testDimension() {
162 1 : const size_t n = 16;
163 1 : const size_t m = 9;
164 1 : const size_t p = 23;
165 1 : Matrix F = MatrixFactory::MakeRandomMatrix(p, n, 0.0, 1.0, Matrix::MATRIX_DENSE);
166 2 : Matrix A = MatrixFactory::MakeRandomMatrix(p, n + 1, 0.0, 1.0, Matrix::MATRIX_DENSE);
167 2 : Matrix B = MatrixFactory::MakeRandomMatrix(n, m, 0.0, 1.0, Matrix::MATRIX_DENSE);
168 1 : MatrixOperator *T1 = new MatrixOperator(A);
169 1 : MatrixOperator *T2 = new MatrixOperator(B);
170 : OpComposition *op;
171 1 : _ASSERT_EXCEPTION(op = new OpComposition(*T1, *T2), std::invalid_argument);
172 :
173 1 : T1 = new MatrixOperator(F);
174 1 : T2 = new MatrixOperator(B);
175 :
176 1 : op = new OpComposition(*T1, *T2);
177 1 : _ASSERT_EQ(m, op->dimensionIn().first);
178 1 : _ASSERT_EQ(p, op->dimensionOut().first);
179 :
180 1 : delete T1;
181 1 : delete T2;
182 2 : delete op;
183 :
184 4 : }
185 :
186 :
187 :
|
