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

Generated on Sun Apr 20 12:40:05 2008 for IT++ by Doxygen 1.5.5