1 /* -*- c++ -*- */
2 /*
3 * Copyright 2015 Free Software Foundation, Inc.
4 *
5 * This is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published
7 * by the Free Software Foundation; either version 3, or (at your
8 * option) any later version.
9 *
10 * This software is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this software; see the file COPYING. If not, write to
17 * the Free Software Foundation, Inc., 51 Franklin Street,
18 * Boston, MA 02110-1301, USA.
19 */
20
21 #ifdef HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24
25 #include "fec_mtrx_impl.h"
26 #include <math.h>
27 #include <fstream>
28 #include <iostream>
29 #include <sstream>
30 #include <stdexcept>
31 #include <vector>
32
33 namespace gr {
34 namespace fec {
35 namespace code {
36
37 // For use when casting to a matrix_sptr to provide the proper
38 // callback to free the memory when the reference count goes to
39 // 0. Needed to know how to cast from our matrix to gsl_matrix.
matrix_free(matrix * x)40 void matrix_free(matrix* x) { gsl_matrix_free((gsl_matrix*)x); }
41
42
read_matrix_from_file(const std::string filename)43 matrix_sptr read_matrix_from_file(const std::string filename)
44 {
45 std::ifstream inputFile;
46 unsigned int ncols, nrows;
47
48 // Open the alist file (needs a C-string)
49 inputFile.open(filename.c_str());
50 if (!inputFile) {
51 std::stringstream s;
52 s << "There was a problem opening file '" << filename << "' for reading.";
53 throw std::runtime_error(s.str());
54 }
55
56 // First line of file is matrix size: # columns, # rows
57 inputFile >> ncols >> nrows;
58
59 // Now we can allocate memory for the GSL matrix
60 gsl_matrix* temp_matrix = gsl_matrix_alloc(nrows, ncols);
61 gsl_matrix_set_zero(temp_matrix);
62
63 // The next few lines in the file are not necessary in
64 // constructing the matrix.
65 std::string tempBuffer;
66 unsigned int counter;
67 for (counter = 0; counter < 4; counter++) {
68 getline(inputFile, tempBuffer);
69 }
70
71 // These next lines list the indices for where the 1s are in
72 // each column.
73 unsigned int column_count = 0;
74 std::string row_number;
75
76 while (column_count < ncols) {
77 getline(inputFile, tempBuffer);
78 std::stringstream ss(tempBuffer);
79
80 while (ss >> row_number) {
81 int row_i = atoi(row_number.c_str());
82 // alist files index starting from 1, not 0, so decrement
83 row_i--;
84 // set the corresponding matrix element to 1
85 gsl_matrix_set(temp_matrix, row_i, column_count, 1);
86 }
87 column_count++;
88 }
89
90 // Close the alist file
91 inputFile.close();
92
93 // Stash the pointer
94 matrix_sptr H = matrix_sptr((matrix*)temp_matrix, matrix_free);
95
96 return H;
97 }
98
write_matrix_to_file(const std::string filename,matrix_sptr M)99 void write_matrix_to_file(const std::string filename, matrix_sptr M)
100 {
101 std::ofstream outputfile;
102
103 // Open the output file
104 outputfile.open(filename.c_str());
105 if (!outputfile) {
106 std::stringstream s;
107 s << "There was a problem opening file '" << filename << "' for writing.";
108 throw std::runtime_error(s.str());
109 }
110
111 unsigned int ncols = M->size2;
112 unsigned int nrows = M->size1;
113 std::vector<unsigned int> colweights(ncols, 0);
114 std::vector<unsigned int> rowweights(nrows, 0);
115 std::stringstream colout;
116 std::stringstream rowout;
117
118 for (unsigned int c = 0; c < ncols; c++) {
119 for (unsigned int r = 0; r < nrows; r++) {
120 double x = gsl_matrix_get((gsl_matrix*)(M.get()), r, c);
121 if (x == 1) {
122 colout << (r + 1) << " ";
123 colweights[c]++;
124 }
125 }
126 colout << std::endl;
127 }
128
129 for (unsigned int r = 0; r < nrows; r++) {
130 for (unsigned int c = 0; c < ncols; c++) {
131 double x = gsl_matrix_get((gsl_matrix*)(M.get()), r, c);
132 if (x == 1) {
133 rowout << (c + 1) << " ";
134 rowweights[r]++;
135 }
136 }
137 rowout << std::endl;
138 }
139
140 outputfile << ncols << " " << nrows << std::endl;
141 outputfile << (*std::max_element(colweights.begin(), colweights.end())) << " "
142 << (*std::max_element(rowweights.begin(), rowweights.end())) << std::endl;
143
144 std::vector<unsigned int>::iterator itr;
145 for (itr = colweights.begin(); itr != colweights.end(); itr++) {
146 outputfile << (*itr) << " ";
147 }
148 outputfile << std::endl;
149
150 for (itr = rowweights.begin(); itr != rowweights.end(); itr++) {
151 outputfile << (*itr) << " ";
152 }
153 outputfile << std::endl;
154
155 outputfile << colout.str() << rowout.str();
156
157 // Close the alist file
158 outputfile.close();
159 }
160
generate_G_transpose(matrix_sptr H_obj)161 matrix_sptr generate_G_transpose(matrix_sptr H_obj)
162 {
163 unsigned int k = H_obj->size1;
164 unsigned int n = H_obj->size2;
165 unsigned int row_index, col_index;
166
167 gsl_matrix* G_transp = gsl_matrix_alloc(n, k);
168
169 // Grab P' matrix (P' denotes P transposed)
170 gsl_matrix* P_transpose = gsl_matrix_alloc(n - k, k);
171 for (row_index = 0; row_index < n - k; row_index++) {
172 for (col_index = 0; col_index < k; col_index++) {
173 int value = gsl_matrix_get((gsl_matrix*)(H_obj.get()), row_index, col_index);
174 gsl_matrix_set(P_transpose, row_index, col_index, value);
175 }
176 }
177
178 // Set G transpose matrix (used for encoding)
179 gsl_matrix_set_zero(G_transp);
180 for (row_index = 0; row_index < k; row_index++) {
181 col_index = row_index;
182 gsl_matrix_set(G_transp, row_index, col_index, 1);
183 }
184 for (row_index = k; row_index < n; row_index++) {
185 for (col_index = 0; col_index < k; col_index++) {
186 int value = gsl_matrix_get(P_transpose, row_index - k, col_index);
187 gsl_matrix_set(G_transp, row_index, col_index, value);
188 }
189 }
190
191 // Stash the pointer
192 matrix_sptr G = matrix_sptr((matrix*)G_transp, matrix_free);
193
194 // Free memory
195 gsl_matrix_free(P_transpose);
196
197 return G;
198 }
199
generate_G(matrix_sptr H_obj)200 matrix_sptr generate_G(matrix_sptr H_obj)
201 {
202 matrix_sptr G_trans = generate_G_transpose(H_obj);
203
204 unsigned int k = H_obj->size1;
205 unsigned int n = H_obj->size2;
206 gsl_matrix* G = gsl_matrix_alloc(k, n);
207
208 gsl_matrix_transpose_memcpy(G, (gsl_matrix*)(G_trans.get()));
209
210 matrix_sptr Gret = matrix_sptr((matrix*)G, matrix_free);
211 return Gret;
212 }
213
generate_H(matrix_sptr G_obj)214 matrix_sptr generate_H(matrix_sptr G_obj)
215 {
216 unsigned int row_index, col_index;
217
218 unsigned int n = G_obj->size2;
219 unsigned int k = G_obj->size1;
220
221 gsl_matrix* G_ptr = (gsl_matrix*)(G_obj.get());
222 gsl_matrix* H_ptr = gsl_matrix_alloc(n - k, n);
223
224 // Grab P matrix
225 gsl_matrix* P = gsl_matrix_alloc(k, n - k);
226 for (row_index = 0; row_index < k; row_index++) {
227 for (col_index = 0; col_index < n - k; col_index++) {
228 int value = gsl_matrix_get(G_ptr, row_index, col_index + k);
229 gsl_matrix_set(P, row_index, col_index, value);
230 }
231 }
232
233 // Calculate P transpose
234 gsl_matrix* P_transpose = gsl_matrix_alloc(n - k, k);
235 gsl_matrix_transpose_memcpy(P_transpose, P);
236
237 // Set H matrix. H = [-P' I] but since we are doing mod 2,
238 // -P = P, so H = [P' I]
239 gsl_matrix_set_zero(H_ptr);
240 for (row_index = 0; row_index < n - k; row_index++) {
241 for (col_index = 0; col_index < k; col_index++) {
242 int value = gsl_matrix_get(P_transpose, row_index, col_index);
243 gsl_matrix_set(H_ptr, row_index, col_index, value);
244 }
245 }
246
247 for (row_index = 0; row_index < n - k; row_index++) {
248 col_index = row_index + k;
249 gsl_matrix_set(H_ptr, row_index, col_index, 1);
250 }
251
252 // Free memory
253 gsl_matrix_free(P);
254 gsl_matrix_free(P_transpose);
255
256 matrix_sptr H = matrix_sptr((matrix*)H_ptr, matrix_free);
257 return H;
258 }
259
260
print_matrix(const matrix_sptr M,bool numpy)261 void print_matrix(const matrix_sptr M, bool numpy)
262 {
263 if (!numpy) {
264 for (size_t i = 0; i < M->size1; i++) {
265 for (size_t j = 0; j < M->size2; j++) {
266 std::cout << gsl_matrix_get((gsl_matrix*)(M.get()), i, j) << " ";
267 }
268 std::cout << std::endl;
269 }
270 std::cout << std::endl;
271 } else {
272 std::cout << "numpy.matrix([ ";
273 for (size_t i = 0; i < M->size1; i++) {
274 std::cout << "[ ";
275 for (size_t j = 0; j < M->size2; j++) {
276 std::cout << gsl_matrix_get((gsl_matrix*)(M.get()), i, j) << ", ";
277 }
278 std::cout << "], ";
279 }
280 std::cout << "])" << std::endl;
281 }
282 }
283
284
fec_mtrx_impl()285 fec_mtrx_impl::fec_mtrx_impl()
286 {
287 // Assume the convention that parity bits come last in the
288 // codeword
289 d_par_bits_last = true;
290 }
291
H() const292 const gsl_matrix* fec_mtrx_impl::H() const
293 {
294 const gsl_matrix* H_ptr = (gsl_matrix*)(d_H_sptr.get());
295 return H_ptr;
296 }
297
n() const298 unsigned int fec_mtrx_impl::n() const { return d_n; }
299
k() const300 unsigned int fec_mtrx_impl::k() const { return d_k; }
301
add_matrices_mod2(gsl_matrix * result,const gsl_matrix * matrix1,const gsl_matrix * matrix2) const302 void fec_mtrx_impl::add_matrices_mod2(gsl_matrix* result,
303 const gsl_matrix* matrix1,
304 const gsl_matrix* matrix2) const
305 {
306 // This function returns ((matrix1 + matrix2) % 2).
307
308 // Verify that matrix sizes are appropriate
309 unsigned int matrix1_rows = (*matrix1).size1;
310 unsigned int matrix1_cols = (*matrix1).size2;
311 unsigned int matrix2_rows = (*matrix2).size1;
312 unsigned int matrix2_cols = (*matrix2).size2;
313
314 if (matrix1_rows != matrix2_rows) {
315 std::cout << "Error in add_matrices_mod2. Matrices do"
316 << " not have the same number of rows.\n";
317 exit(1);
318 }
319 if (matrix1_cols != matrix2_cols) {
320 std::cout << "Error in add_matrices_mod2. Matrices do"
321 << " not have the same number of columns.\n";
322 exit(1);
323 }
324
325 // Copy matrix1 into result
326 gsl_matrix_memcpy(result, matrix1);
327
328 // Do subtraction. This is not mod 2 yet.
329 gsl_matrix_add(result, matrix2);
330
331 // Take care of mod 2 manually
332 unsigned int row_index, col_index;
333 for (row_index = 0; row_index < matrix1_rows; row_index++) {
334 for (col_index = 0; col_index < matrix1_cols; col_index++) {
335 int value = gsl_matrix_get(result, row_index, col_index);
336 int new_value = abs(value) % 2;
337 gsl_matrix_set(result, row_index, col_index, new_value);
338 }
339 }
340 }
341
mult_matrices_mod2(gsl_matrix * result,const gsl_matrix * matrix1,const gsl_matrix * matrix2) const342 void fec_mtrx_impl::mult_matrices_mod2(gsl_matrix* result,
343 const gsl_matrix* matrix1,
344 const gsl_matrix* matrix2) const
345 {
346 // Verify that matrix sizes are appropriate
347 unsigned int a = (*matrix1).size1; // # of rows
348 unsigned int b = (*matrix1).size2; // # of columns
349 unsigned int c = (*matrix2).size1; // # of rows
350 unsigned int d = (*matrix2).size2; // # of columns
351 if (b != c) {
352 std::cout << "Error in "
353 << "fec_mtrx_impl::mult_matrices_mod2."
354 << " Matrix dimensions do not allow for matrix "
355 << "multiplication operation:\nmatrix1 is " << a << " x " << b
356 << ", and matrix2 is " << c << " x " << d << ".\n";
357 exit(1);
358 }
359
360 // Perform matrix multiplication. This is not mod 2.
361 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, matrix1, matrix2, 0.0, result);
362
363 // Take care of mod 2 manually.
364 unsigned int row_index, col_index;
365 unsigned int rows = (*result).size1, cols = (*result).size2;
366 for (row_index = 0; row_index < rows; row_index++) {
367 for (col_index = 0; col_index < cols; col_index++) {
368 int value = gsl_matrix_get(result, row_index, col_index);
369 int new_value = value % 2;
370 gsl_matrix_set(result, row_index, col_index, new_value);
371 }
372 }
373 }
374
calc_inverse_mod2(const gsl_matrix * original_matrix) const375 gsl_matrix* fec_mtrx_impl::calc_inverse_mod2(const gsl_matrix* original_matrix) const
376 {
377
378 // Let n represent the size of the n x n square matrix
379 unsigned int n = (*original_matrix).size1;
380 unsigned int row_index, col_index;
381
382 // Make a copy of the original matrix, call it matrix_altered.
383 // This matrix will be modified by the GSL functions.
384 gsl_matrix* matrix_altered = gsl_matrix_alloc(n, n);
385 gsl_matrix_memcpy(matrix_altered, original_matrix);
386
387 // In order to find the inverse, GSL must perform a LU
388 // decomposition first.
389 gsl_permutation* permutation = gsl_permutation_alloc(n);
390 int signum;
391 gsl_linalg_LU_decomp(matrix_altered, permutation, &signum);
392
393 // Allocate memory to store the matrix inverse
394 gsl_matrix* matrix_inverse = gsl_matrix_alloc(n, n);
395
396 // Find matrix inverse. This is not mod2.
397 int status = gsl_linalg_LU_invert(matrix_altered, permutation, matrix_inverse);
398
399 if (status) {
400 // Inverse not found by GSL functions.
401 throw "Error in calc_inverse_mod2(): inverse not found.\n";
402 }
403
404 // Find determinant
405 float determinant = gsl_linalg_LU_det(matrix_altered, signum);
406
407 // Multiply the matrix inverse by the determinant.
408 gsl_matrix_scale(matrix_inverse, determinant);
409
410 // Take mod 2 of each element in the matrix.
411 for (row_index = 0; row_index < n; row_index++) {
412 for (col_index = 0; col_index < n; col_index++) {
413
414 float value = gsl_matrix_get(matrix_inverse, row_index, col_index);
415
416 // take care of mod 2
417 int value_cast_as_int = static_cast<int>(value);
418 int temp_value = abs(fmod(value_cast_as_int, 2));
419
420 gsl_matrix_set(matrix_inverse, row_index, col_index, temp_value);
421 }
422 }
423
424 int max_value = gsl_matrix_max(matrix_inverse);
425 if (!max_value) {
426 throw "Error in calc_inverse_mod2(): The matrix inverse found is all zeros.\n";
427 }
428
429 // Verify that the inverse was found by taking matrix
430 // product of original_matrix and the inverse, which should
431 // equal the identity matrix.
432 gsl_matrix* test = gsl_matrix_alloc(n, n);
433 mult_matrices_mod2(test, original_matrix, matrix_inverse);
434
435 gsl_matrix* identity = gsl_matrix_alloc(n, n);
436 gsl_matrix_set_identity(identity);
437 // int test_if_equal = gsl_matrix_equal(identity,test);
438 gsl_matrix_sub(identity, test); // should be null set if equal
439
440 double test_if_not_equal = gsl_matrix_max(identity);
441
442 if (test_if_not_equal > 0) {
443 throw "Error in calc_inverse_mod2(): The matrix inverse found is not valid.\n";
444 }
445
446 return matrix_inverse;
447 }
448
parity_bits_come_last() const449 bool fec_mtrx_impl::parity_bits_come_last() const { return d_par_bits_last; }
450
~fec_mtrx_impl()451 fec_mtrx_impl::~fec_mtrx_impl() {}
452
453 } /* namespace code */
454 } /* namespace fec */
455 } /* namespace gr */
456