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