--- /dev/null
+#ifndef booleanmatrix_h_AechohNgakoghahV
+#define booleanmatrix_h_AechohNgakoghahV
+
+#include <stdint.h>
+#include <cassert>
+#include <cstdlib>
+#include <algorithm>
+#include <vector>
+#include <list>
+#include <set>
+#include <ostream>
+#include <iterator>
+
+#include "booleanvector.h"
+
+namespace libstick {
+
+
+template<class T>
+std::ostream& operator<<(std::ostream &, const std::vector<T> &);
+
+
+/** The base class of boolean_colmatrix and boolean_colrowmatrix which implements
+ * the common logic of both. */
+template<class IT, class D>
+class boolean_colmatrix_base {
+
+ public:
+ typedef IT index_type;
+ typedef sorted_boolean_vector<IT> column_type;
+ typedef D derived;
+
+ protected:
+ /** Create a matrix with 'width' columns, initalized with zero entries. */
+ boolean_colmatrix_base(size_t width) :
+ cols(width) {
+ }
+
+ public:
+ /** A casting constructor for any colmatrix with same entries type. */
+ template<class D2>
+ boolean_colmatrix_base(const boolean_colmatrix_base<IT, D2> &mat) :
+ cols(mat.get_columns()) {
+ }
+
+ /** Get width of the matrix. */
+ size_t width() const {
+ return cols.size();
+ }
+
+ /** Get height of the matrix, i.e. maximum row-index + 1 among all
+ * columns. */
+ size_t height() const {
+ IT h = 0;
+ for (unsigned c=0; c < width(); ++c)
+ if (!cols[c].isempty())
+ h = std::max(h, cols[c].back());
+ return h+1;
+ }
+
+ /** Get the matrix entry at row 'r' and column 'c'. */
+ bool get(index_type r, index_type c) const {
+ assert(c < width());
+ return get_column(c).get(r);
+ }
+
+ /** Set the matrix entry at row 'r' and column 'c'. */
+ void set(index_type r, index_type c, bool value) {
+ get_derived()->_set(r, c, value);
+ }
+
+ /** For each of the 'count'-many (row-index, column-pair) pair in
+ * 'indices', set the specific value. */
+ void set_all(index_type indices[][2], size_t count, bool value) {
+ for (unsigned i=0; i < count; ++i)
+ set(indices[i][0], indices[i][1], value);
+ }
+
+ /** Get the c-th column. */
+ const column_type& get_column(index_type c) const {
+ assert(c < width());
+ return cols[c];
+ }
+
+ /** Get all columns */
+ const std::vector<column_type>& get_columns() const {
+ return cols;
+ }
+
+ /** Add the column-vector 'col' to the c-th column. Note that 'col'
+ * actually contains the list of row-indices that are 1. */
+ void add_column(index_type c, const column_type &col) {
+ assert(c < width());
+
+ // Flip all entries that are set in 'col'.
+ for (typename column_type::indexarray::const_iterator it = col.get_ones().begin();
+ it != col.get_ones().end(); ++it)
+ set(*it, c, !get(*it, c));
+ }
+
+ /** Two matrices are equal iff they have the same entries */
+ template<class D2>
+ bool operator==(const boolean_colmatrix_base<IT, D2> &m) const {
+ return cols == m.get_columns();
+ }
+
+ /** Two matrices are equal iff they have the same entries */
+ template<class D2>
+ bool operator!=(const boolean_colmatrix_base<IT, D2> &m) const {
+ return !(*this == m);
+ }
+
+ /** Print index pairs of the 1s as { {r, c}, {r, c}, {r, c} } */
+ void print_indexpairs(std::ostream &os) const {
+ bool first=true;
+
+ os << "{";
+ for (unsigned c=0; c < width(); ++c) {
+ const column_type &col = get_column(c);
+ const typename column_type::indexarray &ones = col.get_ones();
+ typename column_type::indexarray::iterator it = ones.begin();
+
+ for (typename column_type::indexarray::iterator it = ones.begin();
+ it != ones.end(); ++it) {
+ if (first)
+ first = false;
+ else
+ os << ",";
+ os << " {" << *it << ", " << c << "}";
+ }
+ }
+ os << " }";
+ }
+
+ protected:
+ /** Set the matrix entry at row 'r' and column 'c'. */
+ void _set(index_type r, index_type c, bool value) {
+ assert(c < width());
+ cols.at(c).set(r, value);
+ }
+
+ protected:
+ /** The matrix is the set of columns. */
+ std::vector<column_type> cols;
+
+ /** A cast to the derived type */
+ derived* get_derived() {
+ return static_cast<derived*>(this);
+ }
+};
+
+
+/** This is boolean matrix that is a std::vector of column-vectors and in each
+ * column-vector we save the row-indices of the 1s. It is designed to handle a
+ * few 1s and column-wise operations well. */
+template<class IT=uint32_t>
+class boolean_colmatrix : public boolean_colmatrix_base<IT, boolean_colmatrix<IT> > {
+
+ public:
+ typedef IT index_type;
+ typedef boolean_colmatrix_base<IT, boolean_colmatrix<IT> > base;
+ typedef typename base::column_type column_type;
+
+ /** Create a matrix with 'width' columns, initalized with zero entries. */
+ boolean_colmatrix(size_t columns) :
+ base(columns) {
+ }
+
+ /** Set the matrix entry at row 'r' and column 'c'. */
+ void _set(index_type r, index_type c, bool value) {
+ base::_set(r, c, value);
+ }
+
+ /** A faster implementation of boolean_colmatrix_base::add_column(). */
+ void add_column(index_type c, const column_type &col) {
+ assert(c < base::width());
+ base::cols[c].add(col);
+ }
+
+ private:
+};
+
+
+/** The base class of boolean_rowmatrix and boolean_colrowmatrix which implements
+ * the common logic of both. */
+template<class IT, class D>
+class boolean_rowmatrix_base {
+
+ public:
+ typedef IT index_type;
+ typedef sorted_boolean_vector<IT> row_type;
+ typedef D derived;
+
+ protected:
+ /** Create a matrix with 'height' rows, initalized with zero entries.
+ * */
+ boolean_rowmatrix_base(size_t height) :
+ rows(height) {
+ }
+
+ public:
+ /** Get height of the matrix. */
+ size_t height() const {
+ return rows.size();
+ }
+
+ /** Get the matrix entry at row 'r' and column 'c'. */
+ bool get(index_type r, index_type c) const {
+ assert(r < height());
+ return get_row(r).get(c);
+ }
+
+ /** Set the matrix entry at row 'r' and column 'c'. */
+ void set(index_type r, index_type c, bool value) {
+ get_derived()->_set(r, c, value);
+ }
+
+ /** For each of the 'count'-many (row-index, column-pair) pair in
+ * 'indices', set the specific value. */
+ void set_all(index_type indices[][2], size_t count, bool value) {
+ for (unsigned i=0; i < count; ++i)
+ set(indices[i][0], indices[i][1], value);
+ }
+
+ /** Get the r-th row. */
+ const row_type& get_row(index_type r) const {
+ assert(r < height());
+ return rows[r];
+ }
+
+ /** Add the row-vector 'row' to the r-th row. Note that 'row'
+ * actually contains the list of column-indices that are 1. */
+ void add_row(index_type r, const row_type &row) {
+ assert(r < height());
+
+ // Flip all entries that are set in 'row'.
+ for (typename row_type::indexarray::const_iterator it = row.get_ones().begin();
+ it != row.get_ones().end(); ++it)
+ set(r, *it, !get(r, *it));
+ }
+
+ /** Two matrices are equal iff they have the same entries */
+ bool operator==(const boolean_rowmatrix_base<IT, D> &m) const {
+ return rows == m.rows;
+ }
+
+ /** Two matrices are equal iff they have the same entries */
+ bool operator!=(const boolean_rowmatrix_base<IT, D> &m) const {
+ return !(*this == m);
+ }
+
+ protected:
+ /** Set the matrix entry at row 'r' and column 'c'. */
+ void _set(index_type r, index_type c, bool value) {
+ assert(r < height());
+ rows.at(r).set(c, value);
+ }
+
+ protected:
+ derived* get_derived() {
+ return static_cast<derived*>(this);
+ }
+
+ /** The matrix is the set of columns. */
+ std::vector<row_type> rows;
+};
+
+
+/** This is boolean matrix that is a std::vector of row-vectors and in each
+ * row-vector we save the column-indices of the 1s. It is designed to handle a
+ * few 1s and row-wise operations well. */
+template<class IT=uint32_t>
+class boolean_rowmatrix : public boolean_rowmatrix_base<IT, boolean_rowmatrix<IT> > {
+
+ public:
+ typedef IT index_type;
+ typedef boolean_rowmatrix_base<IT, boolean_rowmatrix<IT> > base;
+
+ /** Create a matrix with 'height' rows, initalized with zero entries. */
+ boolean_rowmatrix(size_t height) :
+ base(height) {
+ }
+
+ /** Set the matrix entry at row 'r' and column 'c'. */
+ void _set(index_type r, index_type c, bool value) {
+ base::_set(r, c, value);
+ }
+};
+
+
+/** This is a boolean matrix that supports. It is designed to handle a
+ * few 1s and row- and column-wise operations well. */
+template<class IT=uint32_t>
+class boolean_colrowmatrix : public boolean_colmatrix_base<IT, boolean_colrowmatrix<IT> >,
+ public boolean_rowmatrix_base<IT, boolean_colrowmatrix<IT> > {
+
+ public:
+ typedef boolean_colmatrix_base<IT, boolean_colrowmatrix<IT> > colbase;
+ typedef boolean_rowmatrix_base<IT, boolean_colrowmatrix<IT> > rowbase;
+
+ public:
+ typedef IT index_type;
+
+ /** Create a size x size matrix, initialized with zeros. */
+ boolean_colrowmatrix(size_t size) :
+ colbase(size),
+ rowbase(size) {
+ }
+
+ /** Casting a colmatrix into a colrow matrix */
+ template<class D>
+ boolean_colrowmatrix(const boolean_colmatrix_base<IT, D>& mat) :
+ colbase(std::max(mat.width(), mat.height())),
+ rowbase(std::max(mat.width(), mat.height())) {
+ for (unsigned c=0; c < mat.width(); ++c) {
+ const typename colbase::column_type &col = mat.get_column(c);
+ for (unsigned i=0; i < col.size(); ++i)
+ set(col.get_ones()[i], c, true);
+ }
+ for (unsigned r=0; r < size(); ++r) {
+ const typename rowbase::row_type &row = rowbase::get_row(r);
+ for (unsigned i=0; i < row.size(); ++i)
+ assert(get(r, row.get_ones()[i]) == true);
+ }
+ }
+
+ /** Override implementation. */
+ void _set(index_type r, index_type c, bool value) {
+ colbase::_set(r, c, value);
+ rowbase::_set(r, c, value);
+ }
+
+ public:
+ /** Get height resp. width of the matrix. */
+ size_t size() const {
+ assert(colbase::width() == rowbase::height());
+ return colbase::width();
+ }
+
+ /** Set the matrix entry at row 'r' and column 'c'. */
+ void set(index_type r, index_type c, bool value) {
+ //Calling set() for one base suffices, static polymorphism
+ //calls _set() anyhow.
+ //rowbase::set(r, c, value);
+ colbase::set(r, c, value);
+ }
+
+ /** For each of the 'count'-many (row-index, column-pair) pair in
+ * 'indices', set the specific value. */
+ void set_all(index_type indices[][2], size_t count, bool value) {
+ colbase::set_all(indices, count, value);
+ }
+
+ /** Get the matrix entry at row 'r' and column 'c'. */
+ bool get(index_type r, index_type c) const {
+ assert(colbase::get(r, c) == rowbase::get(r, c));
+ return colbase::get(r, c);
+ }
+
+ /** Two matrices are equal iff they have the same entries */
+ bool operator==(const boolean_colrowmatrix<IT> &m) const {
+ assert(rowbase::operator==(m) == colbase::operator==(m));
+ return colbase::operator==(m);
+ }
+
+ /** Two matrices are equal iff they have the same entries */
+ template<class D>
+ bool operator==(const boolean_colmatrix_base<IT, D> &m) const {
+ return colbase::operator==(m);
+ }
+
+ /** Two matrices are equal iff they have the same entries */
+ bool operator!=(const boolean_colrowmatrix<IT> &m) const {
+ return !(*this == m);
+ }
+
+ /** Multiply with matrix b from the right */
+ template<class D>
+ boolean_colrowmatrix<IT> operator*(const boolean_colmatrix_base<IT, D> &b) const {
+ boolean_colrowmatrix<IT> c(size());
+ multiply_matrix(c, *this, b);
+ return c;
+ }
+};
+
+
+/** Counts the number of common elements in two sorted vectors. Equal to
+ * counting the number of elements given by std::set_intersection. */
+template <class InputIterator1, class InputIterator2>
+size_t count_set_intersection(InputIterator1 first1, InputIterator1 last1, InputIterator2 first2, InputIterator2 last2)
+{
+ size_t count = 0;
+
+ // As long as we did not either end, look for common elements
+ while (first1 != last1 && first2 != last2)
+ {
+ if (*first1 < *first2)
+ ++first1;
+ else if (*first2 < *first1)
+ ++first2;
+ // Element is common to both
+ else {
+ ++count;
+ ++first1;
+ ++first2;
+ }
+ }
+ return count;
+}
+
+/** Multiply a*b and save the product in 'result'. It is assumed that 'result'
+ * is intially empty and has appropriate size. */
+template<class IT, class D1, class D2, class RT>
+void multiply_matrix(RT &result, const boolean_rowmatrix_base<IT, D1> &a, const boolean_colmatrix_base<IT, D2> &b) {
+ assert(a.height() == b.width());
+
+ for (unsigned r=0; r < a.height(); ++r) {
+ const typename boolean_rowmatrix_base<IT, D1>::row_type &row = a.get_row(r);
+ for (unsigned c=0; c < b.width(); ++c) {
+ const typename boolean_colmatrix_base<IT, D2>::column_type &col = b.get_column(c);
+ if (count_set_intersection(
+ row.get_ones().begin(), row.get_ones().end(),
+ col.get_ones().begin(), col.get_ones().end()) % 2 == 1)
+ result.set(r, c, true);
+ }
+ }
+}
+
+template<class MT>
+MT create_unit_matrix(size_t size) {
+ MT mat(size);
+ for (unsigned i=0; i < size; ++i)
+ mat.set(i, i, true);
+ return mat;
+}
+
+template<class IT>
+std::ostream& operator<<(std::ostream &os, const boolean_colrowmatrix<IT> &mat) {
+
+ os << " ";
+ for (unsigned c=0; c < mat.size(); ++c) {
+ const unsigned d = (c % 10);
+ if (d > 0)
+ os << d;
+ else
+ os << " ";
+ }
+ os << std::endl;
+
+ for (unsigned r=0; r < mat.size(); ++r) {
+ const unsigned d = (r % 10);
+ if (d > 0)
+ os << d;
+ else
+ os << " ";
+
+ for (unsigned c=0; c < mat.size(); ++c)
+ os << (mat.get(r,c) ? "X" : ".");
+
+ if (r < mat.size() - 1)
+ os << std::endl;
+ }
+ return os;
+}
+
+template<class IT>
+std::ostream& operator<<(std::ostream &os, boolean_colrowmatrix<IT> &mat) {
+ const boolean_colrowmatrix<IT> &m = mat;
+ return os << m;
+}
+
+template<class IT, class D>
+std::ostream& operator<<(std::ostream &os, const boolean_rowmatrix_base<IT, D> &mat) {
+ for (unsigned r=0; r < mat.height(); ++r) {
+ // Get the sorted r-th row
+ std::list<IT> row(mat.get_row(r).size());
+ copy(mat.get_row(r).begin(), mat.get_row(r).end(), row.begin());
+
+ os << "|";
+ // Print the elements, and remove them
+ for (unsigned c=0; row.size() > 0; ++c) {
+ if (row.front() == c) {
+ os << "X";
+ row.pop_front();
+ } else
+ os << ".";
+ }
+
+ if (r < mat.height() - 1)
+ os << std::endl;
+ }
+ return os;
+}
+
+template<class IT, class D>
+std::ostream& operator<<(std::ostream &os, boolean_rowmatrix_base<IT, D> &mat) {
+ const boolean_rowmatrix_base<IT, D> &m = mat;
+ return os << m;
+}
+
+template<class IT, class D>
+std::ostream& operator<<(std::ostream &os, const boolean_colmatrix_base<IT, D> &mat) {
+ // The sorted columns
+ std::vector<std::list<IT> > cols;
+ // The max height (max. value) of all columns
+ IT height = 0;
+
+ for (unsigned c=0; c < mat.width(); ++c) {
+ // Get the sorted c-th column
+ std::list<IT> col(mat.get_column(c).size());
+ copy(mat.get_column(c).begin(), mat.get_column(c).end(), col.begin());
+ cols.push_back(col);
+
+ os << "-";
+
+ if (!col.isempty())
+ height = std::max(height, col.back());
+ }
+
+ os << std::endl;
+
+ for (unsigned r=0; r <= height; ++r) {
+ // Print the elements, and remove them
+ for (unsigned c=0; c < mat.width(); ++c) {
+ std::list<IT> &col = cols.at(c);
+ // This column is already empty
+ if (col.isempty())
+ os << " ";
+ else if (col.front() == r) {
+ os << "X";
+ col.pop_front();
+ } else
+ os << ".";
+ }
+
+ os << std::endl;
+ }
+ return os;
+}
+
+template<class IT, class D>
+std::ostream& operator<<(std::ostream &os, boolean_colmatrix_base<IT, D> &mat) {
+ const boolean_colmatrix_base<IT, D> &m = mat;
+ return os << m;
+}
+
+template<class T>
+std::ostream& operator<<(std::ostream& os, const std::vector<T> &vec) {
+ os << "[";
+
+ typename std::vector<T>::const_iterator it = vec.begin();
+ while ( it != vec.end()) {
+ os << *it;
+ if (++it != vec.end())
+ os << " ";
+ }
+
+ os << "]";
+ return os;
+}
+
+}
+
+#endif