IT++ Logo
gf2mat.cpp
Go to the documentation of this file.
00001 
00029 #include <itpp/base/gf2mat.h>
00030 #include <itpp/base/specmat.h>
00031 #include <itpp/base/matfunc.h>
00032 #include <itpp/base/converters.h>
00033 #include <iostream>
00034 
00035 namespace itpp
00036 {
00037 
00038 // ====== IMPLEMENTATION OF THE ALIST CLASS ==========
00039 
00040 GF2mat_sparse_alist::GF2mat_sparse_alist(const std::string &fname)
00041     : data_ok(false)
00042 {
00043   read(fname);
00044 }
00045 
00046 void GF2mat_sparse_alist::read(const std::string &fname)
00047 {
00048   std::ifstream file;
00049   std::string line;
00050   std::stringstream ss;
00051 
00052   file.open(fname.c_str());
00053   it_assert(file.is_open(),
00054             "GF2mat_sparse_alist::read(): Could not open file \""
00055             << fname << "\" for reading");
00056 
00057   // parse N and M:
00058   getline(file, line);
00059   ss << line;
00060   ss >> N >> M;
00061   it_assert(!ss.fail(),
00062             "GF2mat_sparse_alist::read(): Wrong alist data (N or M)");
00063   it_assert((N > 0) && (M > 0),
00064             "GF2mat_sparse_alist::read(): Wrong alist data");
00065   ss.seekg(0, std::ios::end);
00066   ss.clear();
00067 
00068   // parse max_num_n and max_num_m
00069   getline(file, line);
00070   ss << line;
00071   ss >> max_num_n >> max_num_m;
00072   it_assert(!ss.fail(),
00073             "GF2mat_sparse_alist::read(): Wrong alist data (max_num_{n,m})");
00074   it_assert((max_num_n >= 0) && (max_num_n <= N) &&
00075             (max_num_m >= 0) && (max_num_m <= M),
00076             "GF2mat_sparse_alist::read(): Wrong alist data");
00077   ss.seekg(0, std::ios::end);
00078   ss.clear();
00079 
00080   // parse weight of each column n
00081   num_nlist.set_size(N);
00082   num_nlist.clear();
00083   getline(file, line);
00084   ss << line;
00085   for (int i = 0; i < N; i++) {
00086     ss >> num_nlist(i);
00087     it_assert(!ss.fail(),
00088               "GF2mat_sparse_alist::read(): Wrong alist data (num_nlist("
00089               << i << "))");
00090     it_assert((num_nlist(i) >= 0) && (num_nlist(i) <= M),
00091               "GF2mat_sparse_alist::read(): Wrong alist data (num_nlist("
00092               << i << "))");
00093   }
00094   ss.seekg(0, std::ios::end);
00095   ss.clear();
00096 
00097   // parse weight of each row m
00098   num_mlist.set_size(M);
00099   num_mlist.clear();
00100   getline(file, line);
00101   ss << line;
00102   for (int i = 0; i < M; i++) {
00103     ss >> num_mlist(i);
00104     it_assert(!ss.fail(),
00105               "GF2mat_sparse_alist::read(): Wrong alist data (num_mlist("
00106               << i << "))");
00107     it_assert((num_mlist(i) >= 0) && (num_mlist(i) <= N),
00108               "GF2mat_sparse_alist::read(): Wrong alist data (num_mlist("
00109               << i << "))");
00110   }
00111   ss.seekg(0, std::ios::end);
00112   ss.clear();
00113 
00114   // parse coordinates in the n direction with non-zero entries
00115   nlist.set_size(N, max_num_n);
00116   nlist.clear();
00117   for (int i = 0; i < N; i++) {
00118     getline(file, line);
00119     ss << line;
00120     for (int j = 0; j < num_nlist(i); j++) {
00121       ss >> nlist(i, j);
00122       it_assert(!ss.fail(),
00123                 "GF2mat_sparse_alist::read(): Wrong alist data (nlist("
00124                 << i << "," << j << "))");
00125       it_assert((nlist(i, j) >= 0) && (nlist(i, j) <= M),
00126                 "GF2mat_sparse_alist::read(): Wrong alist data (nlist("
00127                 << i << "," << j << "))");
00128     }
00129     ss.seekg(0, std::ios::end);
00130     ss.clear();
00131   }
00132 
00133   // parse coordinates in the m direction with non-zero entries
00134   mlist.set_size(M, max_num_m);
00135   mlist.clear();
00136   for (int i = 0; i < M; i++) {
00137     getline(file, line);
00138     ss << line;
00139     for (int j = 0; j < num_mlist(i); j++) {
00140       ss >> mlist(i, j);
00141       it_assert(!ss.fail(),
00142                 "GF2mat_sparse_alist::read(): Wrong alist data (mlist("
00143                 << i << "," << j << "))");
00144       it_assert((mlist(i, j) >= 0) && (mlist(i, j) <= N),
00145                 "GF2mat_sparse_alist::read(): Wrong alist data (mlist("
00146                 << i << "," << j << "))");
00147     }
00148     ss.seekg(0, std::ios::end);
00149     ss.clear();
00150   }
00151 
00152   file.close();
00153   data_ok = true;
00154 }
00155 
00156 void GF2mat_sparse_alist::write(const std::string &fname) const
00157 {
00158   it_assert(data_ok,
00159             "GF2mat_sparse_alist::write(): alist data not ready for writing");
00160 
00161   std::ofstream file(fname.c_str(), std::ofstream::out);
00162   it_assert(file.is_open(),
00163             "GF2mat_sparse_alist::write(): Could not open file \""
00164             << fname << "\" for writing");
00165 
00166   file << N << " " << M << std::endl;
00167   file << max_num_n << " " << max_num_m << std::endl;
00168 
00169   for (int i = 0; i < num_nlist.length() - 1; i++)
00170     file << num_nlist(i) << " ";
00171   file << num_nlist(num_nlist.length() - 1) << std::endl;
00172 
00173   for (int i = 0; i < num_mlist.length() - 1; i++)
00174     file << num_mlist(i) << " ";
00175   file << num_mlist(num_mlist.length() - 1) << std::endl;
00176 
00177   for (int i = 0; i < N; i++) {
00178     for (int j = 0; j < num_nlist(i) - 1; j++)
00179       file << nlist(i, j) << " ";
00180     file << nlist(i, num_nlist(i) - 1) << std::endl;
00181   }
00182 
00183   for (int i = 0; i < M; i++) {
00184     for (int j = 0; j < num_mlist(i) - 1; j++)
00185       file << mlist(i, j) << " ";
00186     file << mlist(i, num_mlist(i) - 1) << std::endl;
00187   }
00188 
00189   file.close();
00190 }
00191 
00192 
00193 GF2mat_sparse GF2mat_sparse_alist::to_sparse(bool transpose) const
00194 {
00195   GF2mat_sparse sbmat(M, N, max_num_m);
00196 
00197   for (int i = 0; i < M; i++) {
00198     for (int j = 0; j < num_mlist(i); j++) {
00199       sbmat.set_new(i, mlist(i, j) - 1, bin(1));
00200     }
00201   }
00202   sbmat.compact();
00203 
00204   if (transpose) {
00205     return sbmat.transpose();
00206   }
00207   else {
00208     return sbmat;
00209   }
00210 }
00211 
00212 
00213 // ----------------------------------------------------------------------
00214 // WARNING: This method is very slow. Sparse_Mat has to be extended with
00215 // some extra functions to improve the performance of this.
00216 // ----------------------------------------------------------------------
00217 void GF2mat_sparse_alist::from_sparse(const GF2mat_sparse &sbmat,
00218                                       bool transpose)
00219 {
00220   if (transpose) {
00221     from_sparse(sbmat.transpose(), false);
00222   }
00223   else {
00224     // check matrix dimension
00225     M = sbmat.rows();
00226     N = sbmat.cols();
00227 
00228     num_mlist.set_size(M);
00229     num_nlist.set_size(N);
00230 
00231     // fill mlist matrix, num_mlist vector and max_num_m
00232     mlist.set_size(0, 0);
00233     int tmp_cols = 0; // initial number of allocated columns
00234     for (int i = 0; i < M; i++) {
00235       ivec temp_row(0);
00236       for (int j = 0; j < N; j++) {
00237         if (sbmat(i, j) == bin(1)) {
00238           temp_row = concat(temp_row, j + 1);
00239         }
00240       }
00241       int trs = temp_row.size();
00242       if (trs > tmp_cols) {
00243         tmp_cols = trs;
00244         mlist.set_size(M, tmp_cols, true);
00245       }
00246       else if (trs < tmp_cols) {
00247         temp_row.set_size(tmp_cols, true);
00248       }
00249       mlist.set_row(i, temp_row);
00250       num_mlist(i) = trs;
00251     }
00252     max_num_m = max(num_mlist);
00253 
00254     // fill nlist matrix, num_nlist vector and max_num_n
00255     nlist.set_size(0, 0);
00256     tmp_cols = 0; // initial number of allocated columns
00257     for (int j = 0; j < N; j++) {
00258       ivec temp_row = sbmat.get_col(j).get_nz_indices() + 1;
00259       int trs = temp_row.size();
00260       if (trs > tmp_cols) {
00261         tmp_cols = trs;
00262         nlist.set_size(N, tmp_cols, true);
00263       }
00264       else if (trs < tmp_cols) {
00265         temp_row.set_size(tmp_cols, true);
00266       }
00267       nlist.set_row(j, temp_row);
00268       num_nlist(j) = trs;
00269     }
00270     max_num_n = max(num_nlist);
00271 
00272     data_ok = true;
00273   }
00274 }
00275 
00276 
00277 // ----------------------------------------------------------------------
00278 // Implementation of a dense GF2 matrix class
00279 // ----------------------------------------------------------------------
00280 
00281 GF2mat::GF2mat(int i, int j): nrows(i), ncols(j),
00282     nwords((j >> shift_divisor) + 1)
00283 {
00284   data.set_size(nrows, nwords);
00285   data.clear();
00286 }
00287 
00288 GF2mat::GF2mat(): nrows(1), ncols(1), nwords(1)
00289 {
00290   data.set_size(nrows, nwords);
00291   data.clear();
00292 }
00293 
00294 GF2mat::GF2mat(const bvec &x, bool is_column)
00295 {
00296   if (is_column) {     // create column vector
00297     nrows = length(x);
00298     ncols = 1;
00299     nwords = 1;
00300     data.set_size(nrows, nwords);
00301     data.clear();
00302     for (int i = 0; i < nrows; i++) {
00303       set(i, 0, x(i));
00304     }
00305   }
00306   else {    // create row vector
00307     nrows = 1;
00308     ncols = length(x);
00309     nwords = (ncols >> shift_divisor) + 1;
00310     data.set_size(nrows, nwords);
00311     data.clear();
00312     for (int i = 0; i < ncols; i++) {
00313       set(0, i, x(i));
00314     }
00315   }
00316 }
00317 
00318 
00319 GF2mat::GF2mat(const bmat &X): nrows(X.rows()), ncols(X.cols())
00320 {
00321   nwords = (ncols >> shift_divisor) + 1;
00322   data.set_size(nrows, nwords);
00323   data.clear();
00324   for (int i = 0; i < nrows; i++) {
00325     for (int j = 0; j < ncols; j++) {
00326       set(i, j, X(i, j));
00327     }
00328   }
00329 }
00330 
00331 
00332 GF2mat::GF2mat(const GF2mat_sparse &X)
00333 {
00334   nrows = X.rows();
00335   ncols = X.cols();
00336   nwords = (ncols >> shift_divisor) + 1;
00337   data.set_size(nrows, nwords);
00338   for (int i = 0; i < nrows; i++) {
00339     for (int j = 0; j < nwords; j++) {
00340       data(i, j) = 0;
00341     }
00342   }
00343 
00344   for (int j = 0; j < ncols; j++) {
00345     for (int i = 0; i < X.get_col(j).nnz(); i++)   {
00346       bin b = X.get_col(j).get_nz_data(i);
00347       set(X.get_col(j).get_nz_index(i), j, b);
00348     }
00349   }
00350 }
00351 
00352 GF2mat::GF2mat(const GF2mat_sparse &X, int m1, int n1, int m2, int n2)
00353 {
00354   it_assert(X.rows() > m2, "GF2mat(): indexes out of range");
00355   it_assert(X.cols() > n2, "GF2mat(): indexes out of range");
00356   it_assert(m1 >= 0 && n1 >= 0 && m2 >= m1 && n2 >= n1,
00357             "GF2mat::GF2mat(): indexes out of range");
00358 
00359   nrows = m2 - m1 + 1;
00360   ncols = n2 - n1 + 1;
00361   nwords = (ncols >> shift_divisor) + 1;
00362   data.set_size(nrows, nwords);
00363 
00364   for (int i = 0; i < nrows; i++) {
00365     for (int j = 0; j < nwords; j++) {
00366       data(i, j) = 0;
00367     }
00368   }
00369 
00370   for (int i = 0; i < nrows; i++) {
00371     for (int j = 0; j < ncols; j++)   {
00372       bin b = X(i + m1, j + n1);
00373       set(i, j, b);
00374     }
00375   }
00376 }
00377 
00378 
00379 GF2mat::GF2mat(const GF2mat_sparse &X, const ivec &columns)
00380 {
00381   it_assert(X.cols() > max(columns),
00382             "GF2mat::GF2mat(): index out of range");
00383   it_assert(min(columns) >= 0,
00384             "GF2mat::GF2mat(): column index must be positive");
00385 
00386   nrows = X.rows();
00387   ncols = length(columns);
00388   nwords = (ncols >> shift_divisor) + 1;
00389   data.set_size(nrows, nwords);
00390 
00391   for (int i = 0; i < nrows; i++) {
00392     for (int j = 0; j < nwords; j++) {
00393       data(i, j) = 0;
00394     }
00395   }
00396 
00397   for (int j = 0; j < ncols; j++) {
00398     for (int i = 0; i < X.get_col(columns(j)).nnz(); i++)   {
00399       bin b = X.get_col(columns(j)).get_nz_data(i);
00400       set(X.get_col(columns(j)).get_nz_index(i), j, b);
00401     }
00402   }
00403 }
00404 
00405 
00406 void GF2mat::set_size(int m, int n, bool copy)
00407 {
00408   nrows = m;
00409   ncols = n;
00410   nwords = (ncols >> shift_divisor) + 1;
00411   data.set_size(nrows, nwords, copy);
00412   if (!copy)
00413     data.clear();
00414 }
00415 
00416 
00417 GF2mat_sparse GF2mat::sparsify() const
00418 {
00419   GF2mat_sparse Z(nrows, ncols);
00420   for (int i = 0; i < nrows; i++) {
00421     for (int j = 0; j < ncols; j++) {
00422       if (get(i, j) == 1) {
00423         Z.set(i, j, 1);
00424       }
00425     }
00426   }
00427 
00428   return Z;
00429 }
00430 
00431 bvec GF2mat::bvecify() const
00432 {
00433   it_assert(nrows == 1 || ncols == 1,
00434             "GF2mat::bvecify() matrix must be a vector");
00435   int n = (nrows == 1 ? ncols : nrows);
00436   bvec result(n);
00437   if (nrows == 1) {
00438     for (int i = 0; i < n; i++) {
00439       result(i) = get(0, i);
00440     }
00441   }
00442   else {
00443     for (int i = 0; i < n; i++) {
00444       result(i) = get(i, 0);
00445     }
00446   }
00447   return result;
00448 }
00449 
00450 
00451 void GF2mat::set_row(int i, bvec x)
00452 {
00453   it_assert(length(x) == ncols,
00454             "GF2mat::set_row(): dimension mismatch");
00455   for (int j = 0; j < ncols; j++) {
00456     set(i, j, x(j));
00457   }
00458 }
00459 
00460 void GF2mat::set_col(int j, bvec x)
00461 {
00462   it_assert(length(x) == nrows,
00463             "GF2mat::set_col(): dimension mismatch");
00464   for (int i = 0; i < nrows; i++) {
00465     set(i, j, x(i));
00466   }
00467 }
00468 
00469 
00470 GF2mat gf2dense_eye(int m)
00471 {
00472   GF2mat Z(m, m);
00473   for (int i = 0; i < m; i++) {
00474     Z.set(i, i, 1);
00475   }
00476   return Z;
00477 }
00478 
00479 GF2mat GF2mat::get_submatrix(int m1, int n1, int m2, int n2) const
00480 {
00481   it_assert(m1 >= 0 && n1 >= 0 && m2 >= m1 && n2 >= n1
00482             && m2 < nrows && n2 < ncols,
00483             "GF2mat::get_submatrix() index out of range");
00484   GF2mat result(m2 - m1 + 1, n2 - n1 + 1);
00485 
00486   for (int i = m1; i <= m2; i++) {
00487     for (int j = n1; j <= n2; j++) {
00488       result.set(i - m1, j - n1, get(i, j));
00489     }
00490   }
00491 
00492   return result;
00493 }
00494 
00495 
00496 GF2mat GF2mat::concatenate_vertical(const GF2mat &X) const
00497 {
00498   it_assert(X.ncols == ncols,
00499             "GF2mat::concatenate_vertical(): dimension mismatch");
00500   it_assert(X.nwords == nwords,
00501             "GF2mat::concatenate_vertical(): dimension mismatch");
00502 
00503   GF2mat result(nrows + X.nrows, ncols);
00504   for (int i = 0; i < nrows; i++) {
00505     for (int j = 0; j < nwords; j++) {
00506       result.data(i, j) = data(i, j);
00507     }
00508   }
00509 
00510   for (int i = 0; i < X.nrows; i++) {
00511     for (int j = 0; j < nwords; j++) {
00512       result.data(i + nrows, j) = X.data(i, j);
00513     }
00514   }
00515 
00516   return result;
00517 }
00518 
00519 GF2mat GF2mat::concatenate_horizontal(const GF2mat &X) const
00520 {
00521   it_assert(X.nrows == nrows,
00522             "GF2mat::concatenate_horizontal(): dimension mismatch");
00523 
00524   GF2mat result(nrows, X.ncols + ncols);
00525   for (int i = 0; i < nrows; i++) {
00526     for (int j = 0; j < ncols; j++) {
00527       result.set(i, j, get(i, j));
00528     }
00529   }
00530 
00531   for (int i = 0; i < nrows; i++) {
00532     for (int j = 0; j < X.ncols; j++) {
00533       result.set(i, j + ncols, X.get(i, j));
00534     }
00535   }
00536 
00537   return result;
00538 }
00539 
00540 bvec GF2mat::get_row(int i) const
00541 {
00542   bvec result(ncols);
00543   for (int j = 0; j < ncols; j++) {
00544     result(j) = get(i, j);
00545   }
00546 
00547   return result;
00548 }
00549 
00550 bvec GF2mat::get_col(int j) const
00551 {
00552   bvec result(nrows);
00553   for (int i = 0; i < nrows; i++) {
00554     result(i) = get(i, j);
00555   }
00556 
00557   return result;
00558 }
00559 
00560 
00561 int GF2mat::T_fact(GF2mat &T, GF2mat &U, ivec &perm) const
00562 {
00563   T = gf2dense_eye(nrows);
00564   U = *this;
00565 
00566   perm = zeros_i(ncols);
00567   for (int i = 0; i < ncols; i++) {
00568     perm(i) = i;
00569   }
00570 
00571   if (nrows > 250) {   // avoid cluttering output ...
00572     it_info_debug("Performing T-factorization of GF(2) matrix...  rows: "
00573                   << nrows << " cols: "  << ncols << " .... " << std::endl);
00574   }
00575   int pdone = 0;
00576   for (int j = 0; j < nrows; j++) {
00577     // Now working on diagonal element j,j
00578     // First try find a row with a 1 in column i
00579     int i1, j1;
00580     for (i1 = j; i1 < nrows; i1++) {
00581       for (j1 = j; j1 < ncols; j1++) {
00582         if (U.get(i1, j1) == 1) { goto found; }
00583       }
00584     }
00585 
00586     return j;
00587 
00588   found:
00589     U.swap_rows(i1, j);
00590     T.swap_rows(i1, j);
00591     U.swap_cols(j1, j);
00592 
00593     int temp = perm(j);
00594     perm(j) = perm(j1);
00595     perm(j1) = temp;
00596 
00597     // now subtract row i from remaining rows
00598     for (int i1 = j + 1; i1 < nrows; i1++)  {
00599       if (U.get(i1, j) == 1) {
00600         int ptemp = floor_i(100.0 * (i1 + j * nrows) / (nrows * nrows));
00601         if (nrows > 250 && ptemp > pdone + 10) {
00602           it_info_debug(ptemp << "% done.");
00603           pdone = ptemp;
00604         }
00605         U.add_rows(i1, j);
00606         T.add_rows(i1, j);
00607       }
00608     }
00609   }
00610   return nrows;
00611 }
00612 
00613 
00614 int GF2mat::T_fact_update_bitflip(GF2mat &T, GF2mat &U,
00615                                   ivec &perm, int rank,
00616                                   int r, int c) const
00617 {
00618   // First, update U (before re-triangulization)
00619   int c0 = c;
00620   for (c = 0; c < ncols; c++) {
00621     if (perm(c) == c0) {
00622       goto foundidx;
00623     }
00624   }
00625   it_error("GF2mat::T_fact_update_bitflip() - internal error");
00626 
00627 foundidx:
00628   for (int i = 0; i < nrows; i++) {
00629     if (T.get(i, r) == 1) {
00630       U.addto_element(i, c, 1);
00631     }
00632   }
00633 
00634   // first move column c to the end
00635   bvec lastcol = U.get_col(c);
00636   int temp_perm = perm(c);
00637   for (int j = c; j < ncols - 1; j++) {
00638     U.set_col(j, U.get_col(j + 1));
00639     perm(j) = perm(j + 1);
00640   }
00641   U.set_col(ncols - 1, lastcol);
00642   perm(ncols - 1) = temp_perm;
00643 
00644   // then, if the matrix is tall, also move row c to the end
00645   if (nrows >= ncols) {
00646     bvec lastrow_U = U.get_row(c);
00647     bvec lastrow_T = T.get_row(c);
00648     for (int i = c; i < nrows - 1; i++) {
00649       U.set_row(i, U.get_row(i + 1));
00650       T.set_row(i, T.get_row(i + 1));
00651     }
00652     U.set_row(nrows - 1, lastrow_U);
00653     T.set_row(nrows - 1, lastrow_T);
00654 
00655     // Do Gaussian elimination on the last row
00656     for (int j = c; j < ncols; j++)  {
00657       if (U.get(nrows - 1, j) == 1) {
00658         U.add_rows(nrows - 1, j);
00659         T.add_rows(nrows - 1, j);
00660       }
00661     }
00662   }
00663 
00664   // Now, continue T-factorization from the point (rank-1,rank-1)
00665   for (int j = rank - 1; j < nrows; j++) {
00666     int i1, j1;
00667     for (i1 = j; i1 < nrows; i1++) {
00668       for (j1 = j; j1 < ncols; j1++) {
00669         if (U.get(i1, j1) == 1) {
00670           goto found;
00671         }
00672       }
00673     }
00674 
00675     return j;
00676 
00677   found:
00678     U.swap_rows(i1, j);
00679     T.swap_rows(i1, j);
00680     U.swap_cols(j1, j);
00681 
00682     int temp = perm(j);
00683     perm(j) = perm(j1);
00684     perm(j1) = temp;
00685 
00686     for (int i1 = j + 1; i1 < nrows; i1++)  {
00687       if (U.get(i1, j) == 1) {
00688         U.add_rows(i1, j);
00689         T.add_rows(i1, j);
00690       }
00691     }
00692   }
00693 
00694   return nrows;
00695 }
00696 
00697 bool GF2mat::T_fact_update_addcol(GF2mat &T, GF2mat &U,
00698                                   ivec &perm, bvec newcol) const
00699 {
00700   int i0 = T.rows();
00701   int j0 = U.cols();
00702   it_assert(j0 > 0, "GF2mat::T_fact_update_addcol(): dimension mismatch");
00703   it_assert(i0 == T.cols(),
00704             "GF2mat::T_fact_update_addcol(): dimension mismatch");
00705   it_assert(i0 == U.rows(),
00706             "GF2mat::T_fact_update_addcol(): dimension mismatch");
00707   it_assert(length(perm) == j0,
00708             "GF2mat::T_fact_update_addcol(): dimension mismatch");
00709   it_assert(U.get(j0 - 1, j0 - 1) == 1,
00710             "GF2mat::T_fact_update_addcol(): dimension mismatch");
00711   // The following test is VERY TIME-CONSUMING
00712   it_assert_debug(U.row_rank() == j0,
00713                   "GF2mat::T_fact_update_addcol(): factorization has incorrect rank");
00714 
00715   bvec z = T * newcol;
00716   GF2mat Utemp = U.concatenate_horizontal(GF2mat(z, 1));
00717 
00718   // start working on position (j0,j0)
00719   int i;
00720   for (i = j0; i < i0; i++) {
00721     if (Utemp.get(i, j0) == 1) {
00722       goto found;
00723     }
00724   }
00725   return (false); // adding the new column would not improve the rank
00726 
00727 found:
00728   perm.set_length(j0 + 1, true);
00729   perm(j0) = j0;
00730 
00731   Utemp.swap_rows(i, j0);
00732   T.swap_rows(i, j0);
00733 
00734   for (int i1 = j0 + 1; i1 < i0; i1++)  {
00735     if (Utemp.get(i1, j0) == 1) {
00736       Utemp.add_rows(i1, j0);
00737       T.add_rows(i1, j0);
00738     }
00739   }
00740 
00741   U = Utemp;
00742   return (true); // the new column was successfully added
00743 }
00744 
00745 
00746 
00747 
00748 GF2mat GF2mat::inverse() const
00749 {
00750   it_assert(nrows == ncols, "GF2mat::inverse(): Matrix must be square");
00751 
00752   // first compute the T-factorization
00753   GF2mat T, U;
00754   ivec perm;
00755   int rank = T_fact(T, U, perm);
00756   it_assert(rank == ncols, "GF2mat::inverse(): Matrix is not full rank");
00757 
00758   // backward substitution
00759   for (int i = ncols - 2; i >= 0; i--) {
00760     for (int j = ncols - 1; j > i; j--) {
00761       if (U.get(i, j) == 1) {
00762         U.add_rows(i, j);
00763         T.add_rows(i, j);
00764       }
00765     }
00766   }
00767   T.permute_rows(perm, 1);
00768   return T;
00769 }
00770 
00771 int GF2mat::row_rank() const
00772 {
00773   GF2mat T, U;
00774   ivec perm;
00775   int rank = T_fact(T, U, perm);
00776   return rank;
00777 }
00778 
00779 bool GF2mat::is_zero() const
00780 {
00781   for (int i = 0; i < nrows; i++) {
00782     for (int j = 0; j < nwords; j++) {
00783       if (data(i, j) != 0) {
00784         return false;
00785       }
00786     }
00787   }
00788   return true;
00789 }
00790 
00791 bool GF2mat::operator==(const GF2mat &X) const
00792 {
00793   if (X.nrows != nrows) { return false; }
00794   if (X.ncols != ncols) { return false; }
00795   it_assert(X.nwords == nwords, "GF2mat::operator==() dimension mismatch");
00796 
00797   for (int i = 0; i < nrows; i++) {
00798     for (int j = 0; j < nwords; j++) {
00799       //      if (X.get(i,j)!=get(i,j)) {
00800       if (X.data(i, j) != data(i, j)) {
00801         return false;
00802       }
00803     }
00804   }
00805   return true;
00806 }
00807 
00808 
00809 void GF2mat::add_rows(int i, int j)
00810 {
00811   it_assert(i >= 0 && i < nrows, "GF2mat::add_rows(): out of range");
00812   it_assert(j >= 0 && j < nrows, "GF2mat::add_rows(): out of range");
00813   for (int k = 0; k < nwords; k++)   {
00814     data(i, k) ^= data(j, k);
00815   }
00816 }
00817 
00818 void GF2mat::swap_rows(int i, int j)
00819 {
00820   it_assert(i >= 0 && i < nrows, "GF2mat::swap_rows(): index out of range");
00821   it_assert(j >= 0 && j < nrows, "GF2mat::swap_rows(): index out of range");
00822   for (int k = 0; k < nwords; k++) {
00823     unsigned char temp = data(i, k);
00824     data(i, k) = data(j, k);
00825     data(j, k) = temp;
00826   }
00827 }
00828 
00829 void GF2mat::swap_cols(int i, int j)
00830 {
00831   it_assert(i >= 0 && i < ncols, "GF2mat::swap_cols(): index out of range");
00832   it_assert(j >= 0 && j < ncols, "GF2mat::swap_cols(): index out of range");
00833   bvec temp = get_col(i);
00834   set_col(i, get_col(j));
00835   set_col(j, temp);
00836 }
00837 
00838 
00839 void GF2mat::operator=(const GF2mat &X)
00840 {
00841   nrows = X.nrows;
00842   ncols = X.ncols;
00843   nwords = X.nwords;
00844   data = X.data;
00845 }
00846 
00847 GF2mat operator*(const GF2mat &X, const GF2mat &Y)
00848 {
00849   it_assert(X.ncols == Y.nrows, "GF2mat::operator*(): dimension mismatch");
00850   it_assert(X.nwords > 0, "Gfmat::operator*(): dimension mismatch");
00851   it_assert(Y.nwords > 0, "Gfmat::operator*(): dimension mismatch");
00852 
00853   /*
00854   // this can be done more efficiently?
00855   GF2mat result(X.nrows,Y.ncols);
00856   for (int i=0; i<X.nrows; i++) {
00857   for (int j=0; j<Y.ncols; j++) {
00858   bin b=0;
00859   for (int k=0; k<X.ncols; k++) {
00860   bin x = X.get(i,k);
00861   bin y = Y.get(k,j);
00862   b ^= (x&y);
00863   }
00864   result.set(i,j,b);
00865   }
00866   }
00867   return result; */
00868 
00869   // is this better?
00870   return mult_trans(X, Y.transpose());
00871 }
00872 
00873 bvec operator*(const GF2mat &X, const bvec &y)
00874 {
00875   it_assert(length(y) == X.ncols, "GF2mat::operator*(): dimension mismatch");
00876   it_assert(X.nwords > 0, "Gfmat::operator*(): dimension mismatch");
00877 
00878   /*
00879   // this can be done more efficiently?
00880   bvec result = zeros_b(X.nrows);
00881   for (int i=0; i<X.nrows; i++) {
00882   // do the nwords-1 data columns first
00883   for (int j=0; j<X.nwords-1; j++) {
00884   int ind = j<<shift_divisor;
00885   unsigned char r=X.data(i,j);
00886   while (r) {
00887   result(i) ^= (r & y(ind));
00888   r >>= 1;
00889   ind++;
00890   }
00891   }
00892   // do the last column separately
00893   for (int j=(X.nwords-1)<<shift_divisor; j<X.ncols; j++) {
00894   result(i) ^= (X.get(i,j) & y(j));
00895   }
00896   }
00897   return result; */
00898 
00899   // is this better?
00900   return (mult_trans(X, GF2mat(y, 0))).bvecify();
00901 }
00902 
00903 GF2mat mult_trans(const GF2mat &X, const GF2mat &Y)
00904 {
00905   it_assert(X.ncols == Y.ncols, "GF2mat::mult_trans(): dimension mismatch");
00906   it_assert(X.nwords > 0, "GF2mat::mult_trans(): dimension mismatch");
00907   it_assert(Y.nwords > 0, "GF2mat::mult_trans(): dimension mismatch");
00908   it_assert(X.nwords == Y.nwords, "GF2mat::mult_trans(): dimension mismatch");
00909 
00910   GF2mat result(X.nrows, Y.nrows);
00911 
00912   for (int i = 0; i < X.nrows; i++) {
00913     for (int j = 0; j < Y.nrows; j++) {
00914       bin b = 0;
00915       int kindx = i;
00916       int kindy = j;
00917       for (int k = 0; k < X.nwords; k++) {
00918         unsigned char r = X.data(kindx) & Y.data(kindy);
00919         /* The following can be speeded up by using a small lookup
00920            table for the number of ones and shift only a few times (or
00921            not at all if table is large) */
00922         while (r) {
00923           b ^= r & 1;
00924           r >>= 1;
00925         };
00926         kindx += X.nrows;
00927         kindy += Y.nrows;
00928       }
00929       result.set(i, j, b);
00930     }
00931   }
00932   return result;
00933 }
00934 
00935 GF2mat GF2mat::transpose() const
00936 {
00937   // CAN BE SPEEDED UP
00938   GF2mat result(ncols, nrows);
00939 
00940   for (int i = 0; i < nrows; i++) {
00941     for (int j = 0; j < ncols; j++) {
00942       result.set(j, i, get(i, j));
00943     }
00944   }
00945   return result;
00946 }
00947 
00948 GF2mat operator+(const GF2mat &X, const GF2mat &Y)
00949 {
00950   it_assert(X.nrows == Y.nrows, "GF2mat::operator+(): dimension mismatch");
00951   it_assert(X.ncols == Y.ncols, "GF2mat::operator+(): dimension mismatch");
00952   it_assert(X.nwords == Y.nwords, "GF2mat::operator+(): dimension mismatch");
00953   GF2mat result(X.nrows, X.ncols);
00954 
00955   for (int i = 0; i < X.nrows; i++) {
00956     for (int j = 0; j < X.nwords; j++) {
00957       result.data(i, j) = X.data(i, j) ^ Y.data(i, j);
00958     }
00959   }
00960 
00961   return result;
00962 }
00963 
00964 void GF2mat::permute_cols(ivec &perm, bool I)
00965 {
00966   it_assert(length(perm) == ncols,
00967             "GF2mat::permute_cols(): dimensions do not match");
00968 
00969   GF2mat temp = (*this);
00970   for (int j = 0; j < ncols; j++) {
00971     if (I == 0) {
00972       set_col(j, temp.get_col(perm(j)));
00973     }
00974     else {
00975       set_col(perm(j), temp.get_col(j));
00976     }
00977   }
00978 }
00979 
00980 void GF2mat::permute_rows(ivec &perm, bool I)
00981 {
00982   it_assert(length(perm) == nrows,
00983             "GF2mat::permute_rows(): dimensions do not match");
00984 
00985   GF2mat temp = (*this);
00986   for (int i = 0; i < nrows; i++) {
00987     if (I == 0) {
00988       for (int j = 0; j < nwords; j++) {
00989         data(i, j) = temp.data(perm(i), j);
00990       }
00991     }
00992     else {
00993       for (int j = 0; j < nwords; j++) {
00994         data(perm(i), j) = temp.data(i, j);
00995       }
00996     }
00997   }
00998 }
00999 
01000 
01001 std::ostream &operator<<(std::ostream &os, const GF2mat &X)
01002 {
01003   int i, j;
01004   os << "---- GF(2) matrix of dimension " << X.nrows << "*" << X.ncols
01005   << " -- Density: " << X.density() << " ----" << std::endl;
01006 
01007   for (i = 0; i < X.nrows; i++) {
01008     os << "      ";
01009     for (j = 0; j < X.ncols; j++) {
01010       os << X.get(i, j) << " ";
01011     }
01012     os << std::endl;
01013   }
01014 
01015   return os;
01016 }
01017 
01018 double GF2mat::density() const
01019 {
01020   int no_of_ones = 0;
01021 
01022   for (int i = 0; i < nrows; i++) {
01023     for (int j = 0; j < ncols; j++) {
01024       no_of_ones += (get(i, j) == 1 ? 1 : 0);
01025     }
01026   }
01027 
01028   return ((double) no_of_ones) / (nrows*ncols);
01029 }
01030 
01031 
01032 it_file &operator<<(it_file &f, const GF2mat &X)
01033 {
01034   // 3 64-bit unsigned words for: nrows, ncols and nwords + rest for char data
01035   uint64_t bytecount = 3 * sizeof(uint64_t)
01036                        + X.nrows * X.nwords * sizeof(char);
01037   f.write_data_header("GF2mat", bytecount);
01038 
01039   f.low_level_write(static_cast<uint64_t>(X.nrows));
01040   f.low_level_write(static_cast<uint64_t>(X.ncols));
01041   f.low_level_write(static_cast<uint64_t>(X.nwords));
01042   for (int i = 0; i < X.nrows; i++) {
01043     for (int j = 0; j < X.nwords; j++) {
01044       f.low_level_write(static_cast<char>(X.data(i, j)));
01045     }
01046   }
01047   return f;
01048 }
01049 
01050 it_ifile &operator>>(it_ifile &f, GF2mat &X)
01051 {
01052   it_file::data_header h;
01053 
01054   f.read_data_header(h);
01055   if (h.type == "GF2mat") {
01056     uint64_t tmp;
01057     f.low_level_read(tmp);
01058     X.nrows = static_cast<int>(tmp);
01059     f.low_level_read(tmp);
01060     X.ncols = static_cast<int>(tmp);
01061     f.low_level_read(tmp);
01062     X.nwords = static_cast<int>(tmp);
01063     X.data.set_size(X.nrows, X.nwords);
01064     for (int i = 0; i < X.nrows; i++) {
01065       for (int j = 0; j < X.nwords; j++) {
01066         char r;
01067         f.low_level_read(r);
01068         X.data(i, j) = static_cast<unsigned char>(r);
01069       }
01070     }
01071   }
01072   else {
01073     it_error("it_ifile &operator>>() - internal error");
01074   }
01075 
01076   return f;
01077 }
01078 
01079 } // namespace itpp
01080 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:29 for IT++ by Doxygen 1.7.4