IT++ Logo

gf2mat.cpp

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

Generated on Thu Apr 24 13:38:58 2008 for IT++ by Doxygen 1.5.5