include/libstick-0.1/booleanmatrix.h

   1 #ifndef booleanmatrix_h_AechohNgakoghahV
   2 #define booleanmatrix_h_AechohNgakoghahV
   3
   4 #include <stdint.h>
   5 #include <cassert>
   6 #include <cstdlib>
   7 #include <algorithm>
   8 #include <vector>
   9 #include <list>
  10 #include <set>
  11 #include <ostream>
  12 #include <iterator>
  13
  14
  15 namespace libstick {
  16
  17
  18 template<class T>
  19 std::ostream& operator<<(std::ostream &, const std::vector<T> &);
  20
  21
  22 /** The base class of boolean_colmatrix and boolean_colrowmatrix which implements
  23  * the common logic of both. */
  24 template<class IT, class D>
  25 class boolean_colmatrix_base {
  26
  27     public:
  28         typedef IT index_type;
  29         typedef std::vector<index_type> column_type;
  30         typedef D derived;
  31
  32     protected:
  33         /** Create a matrix with 'width' columns, initalized with zero entries. */
  34         boolean_colmatrix_base(size_t width) :
  35             cols(width) {
  36         }
  37
  38     public:
  39         /** A casting constructor for any colmatrix with same entries type. */
  40         template<class D2>
  41         boolean_colmatrix_base(const boolean_colmatrix_base<IT, D2> &mat) :
  42             cols(mat.get_columns()) {
  43         }
  44
  45         /** Get width of the matrix. */
  46         size_t width() const {
  47             return cols.size();
  48         }
  49
  50         /** Get height of the matrix, i.e. maximum row-index + 1 among all
  51          * columns. */
  52         size_t height() const {
  53             IT h = 0;
  54             for (unsigned c=0; c < width(); ++c)
  55                 if (cols[c].size() > 0)
  56                     h = std::max(h, cols[c].back());
  57             return h+1;
  58         }
  59
  60         /** Get the matrix entry at row 'r' and column 'c'. */
  61         bool get(index_type r, index_type c) const {
  62             assert(c < width());
  63             const column_type &col = get_column(c);
  64             return binary_search(col.begin(), col.end(), r);
  65         }
  66
  67         /** Set the matrix entry at row 'r' and column 'c'. */
  68         void set(index_type r, index_type c, bool value) {
  69             get_derived()->_set(r, c, value);
  70         }
  71
  72         /** For each of the 'count'-many (row-index, column-pair) pair in
  73          * 'indices', set the specific value. */
  74         void set_all(index_type indices[][2], size_t count, bool value) {
  75             for (unsigned i=0; i < count; ++i)
  76                 set(indices[i][0], indices[i][1], value);
  77         }
  78
  79         /** Get the c-th column. */
  80         const column_type& get_column(index_type c) const {
  81             assert(c < width());
  82             return cols[c];
  83         }
  84
  85         /** Get all columns */
  86         const std::vector<column_type>& get_columns() const {
  87             return cols;
  88         }
  89
  90         /** Add the column-vector 'col' to the c-th column. Note that 'col'
  91          * actually contains the list of row-indices that are 1. */
  92         void add_column(index_type c, const column_type &col) {
  93             assert(c < width());
  94
  95             // Flip all entries that are set in 'col'.
  96             for (typename column_type::const_iterator it = col.begin(); it != col.end(); ++it)
  97                 set(*it, c, !get(*it, c));
  98         }
  99
 100         /** Two matrices are equal iff they have the same entries */
 101         template<class D2>
 102         bool operator==(const boolean_colmatrix_base<IT, D2> &m) const {
 103             return cols == m.get_columns();
 104         }
 105
 106         /** Two matrices are equal iff they have the same entries */
 107         template<class D2>
 108         bool operator!=(const boolean_colmatrix_base<IT, D2> &m) const {
 109             return !(*this == m);
 110         }
 111
 112         /** Print index pairs of the 1s as { {r, c}, {r, c}, {r, c} } */
 113         void print_indexpairs(std::ostream &os) const {
 114             bool first=true;
 115
 116             os << "{";
 117             for (unsigned c=0; c < width(); ++c) {
 118                 const column_type &col = get_column(c);
 119                 for (typename column_type::const_iterator it = col.begin(); it != col.end(); ++it) {
 120                     if (first)
 121                         first = false;
 122                     else
 123                         os << ",";
 124                     os << " {" << *it << ", " << c << "}";
 125                 }
 126             }
 127             os << " }";
 128         }
 129
 130     protected:
 131         /** Set the matrix entry at row 'r' and column 'c'. */
 132         void _set(index_type r, index_type c, bool value) {
 133             assert(c < width());
 134
 135             column_type &col = cols.at(c);
 136             // Let us see where to insert the new element
 137             typename column_type::iterator it = lower_bound(col.begin(), col.end(), r);
 138             bool exists = (it != col.end() && *it == r);
 139             assert(get(r,c) == exists);
 140
 141             // Add 'r' to c-th column
 142             if (value) {
 143                 // r is new, insert it
 144                 if (!exists)
 145                     col.insert(it, r);
 146                 assert(get(r,c));
 147             }
 148             // Remove the element
 149             else {
 150                 if (exists)
 151                     col.erase(it);
 152                 assert(!get(r,c));
 153             }
 154
 155 #ifndef NDEBUG
 156             // C++11 would have is_sorted
 157             for (unsigned i=1; i < col.size(); i++)
 158                 assert(col[i-1] < col[i]);
 159 #endif
 160         }
 161
 162     protected:
 163         /** The matrix is the set of columns. */
 164         std::vector<column_type> cols;
 165
 166         /** A cast to the derived type */
 167         derived* get_derived() {
 168             return static_cast<derived*>(this);
 169         }
 170 };
 171
 172
 173 /** This is boolean matrix that is a std::vector of column-vectors and in each
 174  * column-vector we save the row-indices of the 1s. It is designed to handle a
 175  * few 1s and column-wise operations well. */
 176 template<class IT=uint32_t>
 177 class boolean_colmatrix : public boolean_colmatrix_base<IT, boolean_colmatrix<IT> > {
 178
 179     public:
 180         typedef IT index_type;
 181         typedef boolean_colmatrix_base<IT, boolean_colmatrix<IT> > base;
 182         typedef typename base::column_type column_type;
 183
 184         /** Create a matrix with 'width' columns, initalized with zero entries. */
 185         boolean_colmatrix(size_t columns) :
 186             base(columns) {
 187         }
 188
 189         /** Set the matrix entry at row 'r' and column 'c'. */
 190         void _set(index_type r, index_type c, bool value) {
 191             base::_set(r, c, value);
 192         }
 193
 194         /** A faster implementation of boolean_colmatrix_base::add_column().
 195          * Assumes that 'col' is sorted. */
 196         void add_column(index_type c, const column_type &col) {
 197             assert(c < base::width());
 198
 199 #ifndef NDEBUG
 200             for (unsigned i=1; i < col.size(); ++i)
 201                 assert(col[i-1] < col[i]);
 202 #endif
 203
 204             // The original column
 205             column_type &orig_col = base::cols[c];
 206
 207             // Make target column large enough
 208             const size_t maxsize = orig_col.size() + col.size();
 209             if (tcol.size() < maxsize)
 210                 tcol.resize(maxsize);
 211
 212             // Compute symmetric difference
 213             typename column_type::iterator it = std::set_symmetric_difference(
 214                     orig_col.begin(), orig_col.end(), col.begin(), col.end(), tcol.begin());
 215
 216             // Copy back to the original column
 217             orig_col.resize(it - tcol.begin());
 218             std::copy(tcol.begin(), it, orig_col.begin());
 219         }
 220
 221     private:
 222         /** A temporary container to speed up add_column() */
 223         column_type tcol;
 224 };
 225
 226
 227 /** The base class of boolean_rowmatrix and boolean_colrowmatrix which implements
 228  * the common logic of both. */
 229 template<class IT, class D>
 230 class boolean_rowmatrix_base {
 231
 232     public:
 233         typedef IT index_type;
 234         typedef std::vector<index_type> row_type;
 235         typedef D derived;
 236
 237     protected:
 238         /** Create a matrix with 'height' rows, initalized with zero entries.
 239          * */
 240         boolean_rowmatrix_base(size_t height) :
 241             rows(height) {
 242         }
 243
 244     public:
 245         /** Get height of the matrix. */
 246         size_t height() const {
 247             return rows.size();
 248         }
 249
 250         /** Get the matrix entry at row 'r' and column 'c'. */
 251         bool get(index_type r, index_type c) const {
 252             assert(r < height());
 253             const row_type &row = get_row(r);
 254             return binary_search(row.begin(), row.end(), c);
 255         }
 256
 257         /** Set the matrix entry at row 'r' and column 'c'. */
 258         void set(index_type r, index_type c, bool value) {
 259             get_derived()->_set(r, c, value);
 260         }
 261
 262         /** For each of the 'count'-many (row-index, column-pair) pair in
 263          * 'indices', set the specific value. */
 264         void set_all(index_type indices[][2], size_t count, bool value) {
 265             for (unsigned i=0; i < count; ++i)
 266                 set(indices[i][0], indices[i][1], value);
 267         }
 268
 269         /** Get the r-th row. */
 270         const row_type& get_row(index_type r) const {
 271             assert(r < height());
 272             return rows[r];
 273         }
 274
 275         /** Add the row-vector 'row' to the r-th row. Note that 'row'
 276          * actually contains the list of column-indices that are 1. */
 277         void add_row(index_type r, const row_type &row) {
 278             assert(r < height());
 279
 280             // Flip all entries that are set in 'row'.
 281             for (typename row_type::const_iterator it = row.begin(); it != row.end(); ++it)
 282                 set(r, *it, !get(r, *it));
 283         }
 284
 285         /** Two matrices are equal iff they have the same entries */
 286         bool operator==(const boolean_rowmatrix_base<IT, D> &m) const {
 287             return rows == m.rows;
 288         }
 289
 290         /** Two matrices are equal iff they have the same entries */
 291         bool operator!=(const boolean_rowmatrix_base<IT, D> &m) const {
 292             return !(*this == m);
 293         }
 294
 295     protected:
 296         /** Set the matrix entry at row 'r' and column 'c'. */
 297         void _set(index_type r, index_type c, bool value) {
 298             assert(r < height());
 299
 300             row_type &row = rows.at(r);
 301             // Let us see where to insert/remove the new element
 302             typename row_type::iterator it = lower_bound(row.begin(), row.end(), c);
 303             bool exists = (it != row.end() && *it == c);
 304             assert(get(r,c) == exists);
 305
 306             // Add 'r' to c-th column
 307             if (value) {
 308                 // r is new, insert it
 309                 if (!exists)
 310                     row.insert(it, c);
 311                 assert(get(r,c));
 312             }
 313             // Remove the element
 314             else {
 315                 if (exists)
 316                     row.erase(it);
 317                 assert(!get(r,c));
 318             }
 319
 320 #ifndef NDEBUG
 321             // C++11 would have is_sorted
 322             for (unsigned i=1; i < row.size(); i++)
 323                 assert(row[i-1] < row[i]);
 324 #endif
 325         }
 326
 327     protected:
 328         derived* get_derived() {
 329             return static_cast<derived*>(this);
 330         }
 331
 332         /** The matrix is the set of columns. */
 333         std::vector<row_type> rows;
 334 };
 335
 336
 337 /** This is boolean matrix that is a std::vector of row-vectors and in each
 338  * row-vector we save the column-indices of the 1s. It is designed to handle a
 339  * few 1s and row-wise operations well. */
 340 template<class IT=uint32_t>
 341 class boolean_rowmatrix : public boolean_rowmatrix_base<IT, boolean_rowmatrix<IT> > {
 342
 343     public:
 344         typedef IT index_type;
 345         typedef boolean_rowmatrix_base<IT, boolean_rowmatrix<IT> > base;
 346
 347         /** Create a matrix with 'height' rows, initalized with zero entries. */
 348         boolean_rowmatrix(size_t height) :
 349             base(height) {
 350         }
 351
 352         /** Set the matrix entry at row 'r' and column 'c'. */
 353         void _set(index_type r, index_type c, bool value) {
 354             base::_set(r, c, value);
 355         }
 356 };
 357
 358
 359 /** This is a boolean matrix that supports. It is designed to handle a
 360  * few 1s and row- and column-wise operations well. */
 361 template<class IT=uint32_t>
 362 class boolean_colrowmatrix : public boolean_colmatrix_base<IT, boolean_colrowmatrix<IT> >,
 363                             public boolean_rowmatrix_base<IT, boolean_colrowmatrix<IT> > {
 364
 365     public:
 366         typedef boolean_colmatrix_base<IT, boolean_colrowmatrix<IT> > colbase;
 367         typedef boolean_rowmatrix_base<IT, boolean_colrowmatrix<IT> > rowbase;
 368
 369     public:
 370         typedef IT index_type;
 371
 372         /** Create a size x size matrix, initialized with zeros. */
 373         boolean_colrowmatrix(size_t size) :
 374             colbase(size),
 375             rowbase(size) {
 376         }
 377
 378         /** Casting a colmatrix into a colrow matrix */
 379         template<class D>
 380         boolean_colrowmatrix(const boolean_colmatrix_base<IT, D>& mat) :
 381             colbase(std::max(mat.width(), mat.height())),
 382             rowbase(std::max(mat.width(), mat.height())) {
 383             for (unsigned c=0; c < mat.width(); ++c) {
 384                 const typename colbase::column_type &col = mat.get_column(c);
 385                 for (unsigned i=0; i < col.size(); ++i)
 386                     set(col[i], c, true);
 387             }
 388             for (unsigned r=0; r < size(); ++r) {
 389                 const typename rowbase::row_type &row = rowbase::get_row(r);
 390                 for (unsigned i=0; i < row.size(); ++i)
 391                     assert(get(r, row[i]) == true);
 392             }
 393         }
 394
 395         /** Override implementation. */
 396         void _set(index_type r, index_type c, bool value) {
 397             colbase::_set(r, c, value);
 398             rowbase::_set(r, c, value);
 399         }
 400
 401     public:
 402         /** Get height resp. width of the matrix. */
 403         size_t size() const {
 404             assert(colbase::width() == rowbase::height());
 405             return colbase::width();
 406         }
 407
 408         /** Set the matrix entry at row 'r' and column 'c'. */
 409         void set(index_type r, index_type c, bool value) {
 410             //Calling set() for one base suffices, static polymorphism
 411             //calls _set() anyhow.
 412             //rowbase::set(r, c, value);
 413             colbase::set(r, c, value);
 414         }
 415
 416         /** For each of the 'count'-many (row-index, column-pair) pair in
 417          * 'indices', set the specific value. */
 418         void set_all(index_type indices[][2], size_t count, bool value) {
 419             colbase::set_all(indices, count, value);
 420         }
 421
 422         /** Get the matrix entry at row 'r' and column 'c'. */
 423         bool get(index_type r, index_type c) const {
 424             assert(colbase::get(r, c) == rowbase::get(r, c));
 425             return colbase::get(r, c);
 426         }
 427
 428         /** Two matrices are equal iff they have the same entries */
 429         bool operator==(const boolean_colrowmatrix<IT> &m) const {
 430             assert(rowbase::operator==(m) == colbase::operator==(m));
 431             return colbase::operator==(m);
 432         }
 433
 434         /** Two matrices are equal iff they have the same entries */
 435         template<class D>
 436         bool operator==(const boolean_colmatrix_base<IT, D> &m) const {
 437             return colbase::operator==(m);
 438         }
 439
 440         /** Two matrices are equal iff they have the same entries */
 441         bool operator!=(const boolean_colrowmatrix<IT> &m) const {
 442             return !(*this == m);
 443         }
 444
 445         /** Multiply with matrix b from the right */
 446         template<class D>
 447         boolean_colrowmatrix<IT> operator*(const boolean_colmatrix_base<IT, D> &b) const {
 448             boolean_colrowmatrix<IT> c(size());
 449             multiply_matrix(c, *this, b);
 450             return c;
 451         }
 452 };
 453
 454
 455 /** Counts the number of common elements in two sorted vectors. Equal to
 456  * counting the number of elements given by std::set_intersection. */
 457 template <class InputIterator1, class InputIterator2>
 458 size_t count_set_intersection (InputIterator1 first1, InputIterator1 last1, InputIterator2 first2, InputIterator2 last2)
 459 {
 460     size_t count = 0;
 461
 462     // As long as we did not either end, look for common elements
 463     while (first1 != last1 && first2 != last2)
 464     {
 465         if (*first1 < *first2)
 466             ++first1;
 467         else if (*first2 < *first1)
 468             ++first2;
 469         // Element is common to both
 470         else {
 471             ++count;
 472             ++first1;
 473             ++first2;
 474         }
 475     }
 476     return count;
 477 }
 478
 479 /** Multiply a*b and save the product in 'result'. It is assumed that 'result'
 480  * is intially empty and has appropriate size. */
 481 template<class IT, class D1, class D2, class RT>
 482 void multiply_matrix(RT &result, const boolean_rowmatrix_base<IT, D1> &a, const boolean_colmatrix_base<IT, D2> &b) {
 483     assert(a.height() == b.width());
 484
 485     for (unsigned r=0; r < a.height(); ++r) {
 486         const typename boolean_rowmatrix_base<IT, D1>::row_type &row = a.get_row(r);
 487         for (unsigned c=0; c < b.width(); ++c) {
 488             const typename boolean_colmatrix_base<IT, D2>::column_type &col = b.get_column(c);
 489             if (count_set_intersection(row.begin(), row.end(), col.begin(), col.end()) % 2 == 1)
 490                 result.set(r, c, true);
 491         }
 492     }
 493 }
 494
 495 template<class MT>
 496 MT create_unit_matrix(size_t size) {
 497     MT mat(size);
 498     for (unsigned i=0; i < size; ++i)
 499         mat.set(i, i, true);
 500     return mat;
 501 }
 502
 503 template<class IT>
 504 std::ostream& operator<<(std::ostream &os, const boolean_colrowmatrix<IT> &mat) {
 505
 506     os << " ";
 507     for (unsigned c=0; c < mat.size(); ++c) {
 508         const unsigned d = (c % 10);
 509         if (d > 0)
 510             os << d;
 511         else
 512             os << " ";
 513     }
 514     os << std::endl;
 515
 516     for (unsigned r=0; r < mat.size(); ++r) {
 517         const unsigned d = (r % 10);
 518         if (d > 0)
 519             os << d;
 520         else
 521             os << " ";
 522
 523         for (unsigned c=0; c < mat.size(); ++c)
 524             os << (mat.get(r,c) ? "X" : ".");
 525
 526         if (r < mat.size() - 1)
 527             os << std::endl;
 528     }
 529     return os;
 530 }
 531
 532 template<class IT>
 533 std::ostream& operator<<(std::ostream &os, boolean_colrowmatrix<IT> &mat) {
 534     const boolean_colrowmatrix<IT> &m = mat;
 535     return os << m;
 536 }
 537
 538 template<class IT, class D>
 539 std::ostream& operator<<(std::ostream &os, const boolean_rowmatrix_base<IT, D> &mat) {
 540     for (unsigned r=0; r < mat.height(); ++r) {
 541         // Get the sorted r-th row
 542         std::list<IT> row(mat.get_row(r).size());
 543         copy(mat.get_row(r).begin(), mat.get_row(r).end(), row.begin());
 544
 545         os << "|";
 546         // Print the elements, and remove them
 547         for (unsigned c=0; row.size() > 0; ++c) {
 548             if (row.front() == c) {
 549                 os << "X";
 550                 row.pop_front();
 551             } else
 552                 os << ".";
 553         }
 554
 555         if (r < mat.height() - 1)
 556             os << std::endl;
 557     }
 558     return os;
 559 }
 560
 561 template<class IT, class D>
 562 std::ostream& operator<<(std::ostream &os, boolean_rowmatrix_base<IT, D> &mat) {
 563     const boolean_rowmatrix_base<IT, D> &m = mat;
 564     return os << m;
 565 }
 566
 567 template<class IT, class D>
 568 std::ostream& operator<<(std::ostream &os, const boolean_colmatrix_base<IT, D> &mat) {
 569     // The sorted columns
 570     std::vector<std::list<IT> > cols;
 571     // The max height (max. value) of all columns
 572     IT height = 0;
 573
 574     for (unsigned c=0; c < mat.width(); ++c) {
 575         // Get the sorted c-th column
 576         std::list<IT> col(mat.get_column(c).size());
 577         copy(mat.get_column(c).begin(), mat.get_column(c).end(), col.begin());
 578         cols.push_back(col);
 579
 580         os << "-";
 581
 582         if (col.size() != 0)
 583             height = std::max(height, col.back());
 584     }
 585
 586     os << std::endl;
 587
 588     for (unsigned r=0; r <= height; ++r) {
 589         // Print the elements, and remove them
 590         for (unsigned c=0; c < mat.width(); ++c) {
 591             std::list<IT> &col = cols.at(c);
 592             // This column is already empty
 593             if (col.size() == 0)
 594                 os << " ";
 595             else if (col.front() == r) {
 596                 os << "X";
 597                 col.pop_front();
 598             } else
 599                 os << ".";
 600         }
 601
 602         os << std::endl;
 603     }
 604     return os;
 605 }
 606
 607 template<class IT, class D>
 608 std::ostream& operator<<(std::ostream &os, boolean_colmatrix_base<IT, D> &mat) {
 609     const boolean_colmatrix_base<IT, D> &m = mat;
 610     return os << m;
 611 }
 612
 613 template<class T>
 614 std::ostream& operator<<(std::ostream& os, const std::vector<T> &vec) {
 615     os << "[";
 616
 617     typename std::vector<T>::const_iterator it = vec.begin();
 618     while ( it != vec.end()) {
 619         os << *it;
 620         if (++it != vec.end())
 621             os << " ";
 622     }
 623
 624     os << "]";
 625     return os;
 626 }
 627
 628 }
 629
 630 #endif