From 929c696bfb48aa0f96ce2b70b1b926b5f2b2dede Mon Sep 17 00:00:00 2001 From: Stefan Huber Date: Mon, 11 Nov 2013 13:39:15 +0100 Subject: [PATCH] Add matrix reduction code --- CMakeLists.txt | 2 +- include/libstick-0.1/booleanmatrix.h | 62 +++++++++++- include/libstick-0.1/matrixreduction.h | 51 ++++++++++ include/libstick-0.1/simplicialcomplex.h | 64 ++++++++---- tests/booleanmatrix.h | 4 +- tests/simplicialcomplex.h | 118 ++++++++++++++--------- 6 files changed, 230 insertions(+), 71 deletions(-) create mode 100644 include/libstick-0.1/matrixreduction.h diff --git a/CMakeLists.txt b/CMakeLists.txt index b5116d9..1242e40 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ project (libstick) if(NOT CMAKE_SYSTEM_NAME STREQUAL "Windows") - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -ansi -pedantic") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -std=c++98 -pedantic") # I would like to use -pedantic, but VTK devs are not sufficiently # pedantic set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Werror") diff --git a/include/libstick-0.1/booleanmatrix.h b/include/libstick-0.1/booleanmatrix.h index ab368d2..894892f 100644 --- a/include/libstick-0.1/booleanmatrix.h +++ b/include/libstick-0.1/booleanmatrix.h @@ -55,7 +55,7 @@ class BooleanColMatrix_base { /** 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) { + void addColumn(index_type c, const column_type &col) { assert(c < width()); // Flip all entries that are set in 'col'. @@ -101,6 +101,7 @@ class BooleanColMatrix_base { } #ifndef NDEBUG + // C++11 would have is_sorted for (unsigned i=1; i < col.size(); i++) assert(col[i-1] < col[i]); #endif @@ -180,7 +181,7 @@ class BooleanRowMatrix_base { /** 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) { + void addRow(index_type r, const row_type &row) { assert(r < height()); // Flip all entries that are set in 'row'. @@ -224,6 +225,7 @@ class BooleanRowMatrix_base { } #ifndef NDEBUG + // C++11 would have is_sorted for (unsigned i=1; i < row.size(); i++) assert(row[i-1] < row[i]); #endif @@ -317,8 +319,64 @@ class BooleanColRowMatrix : public BooleanColMatrix_base &m) const { return !(*this == m); } + + /** Multiply with matrix b from the right */ + template + BooleanColRowMatrix operator*(const BooleanColMatrix_base &b) { + BooleanColRowMatrix 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 +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 +void multiply_matrix(RT &result, const BooleanRowMatrix_base &a, const BooleanColMatrix_base &b) { + assert(a.height() == b.width()); + + for (unsigned r=0; r < a.height(); ++r) { + const typename BooleanRowMatrix_base::row_type &row = a.getRow(r); + for (unsigned c=0; c < b.width(); ++c) { + const typename BooleanColMatrix_base::column_type &col = b.getColumn(c); + if (count_set_intersection(row.begin(), row.end(), col.begin(), col.end()) % 2 == 1) + result.set(r, c, true); + } + } +} + +template +MT unitMatrix(size_t size) { + MT mat(size); + for (unsigned i=0; i < size; ++i) + mat.set(i, i, true); + return mat; +} + template std::ostream& operator<<(std::ostream &os, const BooleanColRowMatrix &mat) { for (unsigned r=0; r < mat.size(); ++r) { diff --git a/include/libstick-0.1/matrixreduction.h b/include/libstick-0.1/matrixreduction.h new file mode 100644 index 0000000..60d8397 --- /dev/null +++ b/include/libstick-0.1/matrixreduction.h @@ -0,0 +1,51 @@ +#ifndef matrixreduction_h_aPeinoXeiFeethuz +#define matrixreduction_h_aPeinoXeiFeethuz + +#include + +#include "booleanmatrix.h" + + +namespace libstick { + +/** Takes the boundary matrix b and transforms it inplace into a reduced matrix + * b. The matrix v is required to be a unit matrix having the same dimensions + * as b. It is maintained in order to leave the product b*inverse(v) as an + * invariant. Hence, the resulting b is equal to the product of the boundary + * matrix times v. */ +template +void reduceBoundaryMatrix(BooleanColRowMatrix &b, BooleanColMatrix_base &v) { + assert(b.size() == v.width()); + + // Make every column reduced, i.e., it is a zero-column or contains only one 1. + for (unsigned c=0; c < b.size(); ++c) { + const typename BooleanColRowMatrix::colbase::column_type &col = b.getColumn(c); + if (col.size() == 0) + continue; + + // The column is not a zero-column, take the lowest 1, and reduce the other columns at the row + IT r = col.back(); + assert(b.get(r, c)); + + // Get a copy of the r-th row + typename BooleanColRowMatrix::rowbase::row_type row(b.getRow(r).size()); + copy(b.getRow(r).begin(), b.getRow(r).end(), row.begin()); + assert(row.size() >= 1); + + // Get rid of 1s at that row right of column c. + typename BooleanColRowMatrix::row_type::const_iterator it = lower_bound(row.begin(), row.end(), c); + for(; it != row.end(); ++it) { + if (*it <= c) + continue; + + assert(b.get(r, *it)); + b.addColumn(*it, col); + v.addColumn(*it, v.getColumn(c)); + assert(!b.get(r, *it)); + } + } +} + +} + +#endif diff --git a/include/libstick-0.1/simplicialcomplex.h b/include/libstick-0.1/simplicialcomplex.h index d829721..07080ed 100644 --- a/include/libstick-0.1/simplicialcomplex.h +++ b/include/libstick-0.1/simplicialcomplex.h @@ -6,6 +6,7 @@ #include #include #include +#include #include @@ -17,8 +18,12 @@ namespace libstick { /** A simplicial complex is a std::vector of simplices such that each face is * also part of the complex. Every simplex has dimension at most MAXDIM. The * indices of simplices resp. their faces are of type IT. To each simplex a - * value is assigend, which is of type VT. */ -template + * value is assigend, which is of type VT. When a SimplicialComplex is + * instantiated, a single (-1) dimensional simplex is automatically created. + * Each 0-dimensional simplex automatically has this simplex as its face. + * Consequently, the innner class SimplexOrder gives the extended boundary + * matrix. */ +template class SimplicialComplex { public: @@ -33,34 +38,51 @@ class SimplicialComplex { /** A simplex of the complex. */ struct Simplex { /** Dimension of the simplex. */ - size_t dim; + int dim; /** The indices of the faces of the simplex. */ index_type faces[MAXDIM+1]; /** The value of the simplex. */ value_type value; /** Create a new simplex with dimension 'dim', (dim+1)-faces and - * its value. */ - static Simplex create(size_t dim, index_type* faces, value_type value) { - assert(dim <= MAXDIM); + * its value. If simpley is 0-dimensional, its face is + * automatically set to one (-1)-dimensional simplex. */ + static Simplex create(int dim, index_type* faces, value_type value) { + assert(0 <= dim && dim <= MAXDIM); Simplex s; s.dim = dim; s.value = value; - memcpy(s.faces, faces, face_count_bydim(dim)*sizeof(index_type)); + + if (dim > 0) + memcpy(s.faces, faces, faceCountByDim(dim)*sizeof(index_type)); + else + s.faces[0] = 0; + + return s; + } + + /** Create a (-1)-dimensional simplex. It has the lowest possible value. */ + static Simplex create_minusonedim_simplex() { + Simplex s; + + s.dim = -1; + s.faces[0] = 0; + s.value = std::numeric_limits::has_infinity + ? -std::numeric_limits::infinity() + : std::numeric_limits::min(); return s; } /** Get number of faces. */ - size_t face_count() const { - return face_count_bydim(dim); + size_t faceCount() const { + return faceCountByDim(dim); } /** Get number of faces of a dim-dimensional simplex. */ - static size_t face_count_bydim(size_t dim) { - if (dim == 0) - return 0; + static size_t faceCountByDim(int dim) { + assert(-1 <= dim && dim <= MAXDIM); return dim + 1; } }; @@ -104,7 +126,7 @@ class SimplicialComplex { assert(size() == c.size()); for (unsigned i=0; i < size(); ++i) - for (unsigned f=0; f < getSimplex(i).face_count(); ++f) + for (unsigned f=0; f < getSimplex(i).faceCount(); ++f) if (revorder[getSimplex(i).faces[f]] >= i) return false; @@ -142,7 +164,7 @@ class SimplicialComplex { BoundaryMatrix mat(size()); for (unsigned c=0; c < size(); ++c) - for(unsigned r=0; r < getSimplex(c).face_count(); ++r) + for(unsigned r=0; r < getSimplex(c).faceCount(); ++r) mat.set(revorder[getSimplex(c).faces[r]], c, true); return mat; @@ -171,6 +193,11 @@ class SimplicialComplex { }; public: + SimplicialComplex() { + // Add the one minus-one dimensional simplex + addSimplex(Simplex::create_minusonedim_simplex()); + } + /** Return number of simplices. */ size_t size() const { return simplices.size(); @@ -178,7 +205,7 @@ class SimplicialComplex { /** Add a simplex to the complex. The dimension of the faces must be * dim-1, and they must already be part of the complex. */ - void addSimplex(size_t dim, index_type* faces, value_type value) { + void addSimplex(int dim, index_type* faces, value_type value) { addSimplex(Simplex::create(dim, faces, value)); } @@ -186,7 +213,7 @@ class SimplicialComplex { * dim-1, and they must already be part of the complex. */ void addSimplex(Simplex s) { // Check requirements for faces - for (unsigned i=0; i < s.face_count(); ++i) { + for (unsigned i=0; i < s.faceCount(); ++i) { // Faces are already in complex. assert(s.faces[i] < size()); // Faces have dimension dim-1 @@ -199,11 +226,10 @@ class SimplicialComplex { /** Return true iff for each simplex i with dimension dim it holds that * the faces of i are contained in the complex and have dimension dim-1. */ bool isComplex() const { - typename std::vector::const_iterator it = ++simplices.begin(); for (unsigned i=0; i < size(); ++i) { const Simplex &s = simplices[i]; - for (unsigned f=0; f < s.face_count(); ++f) { + for (unsigned f=0; f < s.faceCount(); ++f) { if (s.faces[f] >= size()) return false; @@ -224,7 +250,7 @@ class SimplicialComplex { typename std::vector::const_iterator it = ++simplices.begin(); for (; it != simplices.end(); ++it) - for (unsigned f=0; f < it->face_count(); ++f) + for (unsigned f=0; f < it->faceCount(); ++f) if (simplices[it->faces[f]].value > it->value) return false; diff --git a/tests/booleanmatrix.h b/tests/booleanmatrix.h index 8554bf3..0cba019 100644 --- a/tests/booleanmatrix.h +++ b/tests/booleanmatrix.h @@ -69,7 +69,7 @@ class BooleanmatrixTestSuite: public Test::Suite { col.push_back(0); // Add column and test for values - mat.add_column(5, col); + mat.addColumn(5, col); for (unsigned c=0; c < size; ++c) { if (c==5) { for (unsigned r=0; r < size; ++r ) @@ -96,7 +96,7 @@ class BooleanmatrixTestSuite: public Test::Suite { row.push_back(3); // Add row and test for values - mat.add_row(5, row); + mat.addRow(5, row); for (unsigned r=0; r < size; ++r) { if (r==5) { for (unsigned c=0; c < size; ++c ) diff --git a/tests/simplicialcomplex.h b/tests/simplicialcomplex.h index fa829e1..afc9fab 100644 --- a/tests/simplicialcomplex.h +++ b/tests/simplicialcomplex.h @@ -5,6 +5,7 @@ #include #include +#include using namespace libstick; @@ -32,6 +33,7 @@ class SimplicialComplexTestSuite: public Test::Suite { TEST_ADD(SimplicialComplexTestSuite::test_isOrderFiltration); TEST_ADD(SimplicialComplexTestSuite::test_isOrderMonotone); TEST_ADD(SimplicialComplexTestSuite::test_boundaryMatrix); + TEST_ADD(SimplicialComplexTestSuite::test_matrixreduction); } protected: @@ -43,43 +45,43 @@ class SimplicialComplexTestSuite: public Test::Suite { const unsigned num = 11; scomplex::Simplex ss[num] = { // dimension, faces, value... - {0, {0, 0, 0}, 0}, - {0, {0, 0, 0}, 1}, - {0, {0, 0, 0}, 2}, - {0, {0, 0, 0}, 3}, - {1, {0, 1, 0}, 4}, - {1, {1, 2, 0}, 5}, - {1, {2, 3, 0}, 6}, - {1, {3, 0, 0}, 7}, - {1, {0, 2, 0}, 8}, - {2, {6, 7, 8}, 9}, - {2, {4, 5, 8}, 10} + {0, {}, 1}, + {0, {}, 2}, + {0, {}, 3}, + {0, {}, 4}, + {1, {1, 2}, 5}, + {1, {2, 3}, 6}, + {1, {3, 4}, 7}, + {1, {4, 1}, 8}, + {1, {1, 3}, 9}, + {2, {7, 8, 9}, 10}, + {2, {5, 6, 9}, 11} }; // This is o1 This is o2 This is o3(b) // (value = index) (values) (values) // - // 0 ----4---- 1 0 ----4---- 1 0 ----4---- 1 - // |\ | |\ | |\ | - // | \ 10 | | \ 10 | | \ 12 | - // | \ | | \ | | \ | - // | \ | | \ | | \ | - // 7 8 5 7 8 11 7 8 11 - // | \ | | \ | | \ | - // | 9 \ | | 9 \ | | 9 \ | - // | \ | | \ | | \ | - // | \| | \| | \| - // 3 ----6---- 2 3 ----6---- 2 3 ----6---- 2 + // 1 ----5---- 2 1 ----5---- 2 1 ----5---- 2 + // |\ | |\ | |\ | + // | \ 11 | | \ 11 | | \ 13 | + // | \ | | \ | | \ | + // | \ | | \ | | \ | + // 8 9 6 8 9 12 8 9 12 + // | \ | | \ | | \ | + // | 10 \ | | 10 \ | | 10 \ | + // | \ | | \ | | \ | + // | \| | \| | \| + // 4 ----7---- 3 4 ----7---- 3 4 ----7---- 3 // Build the complex for (unsigned i=0; i(b.size()); + + bm b_orig = b; + reduceBoundaryMatrix(b, v); + TEST_ASSERT(b == b_orig*v); + + //std::cout << std::endl << "after reduce: " << std::endl; + //std::cout << std::endl << b << std::endl; + //std::cout << std::endl << v << std::endl; + //std::cout << std::endl << (b_orig*v) << std::endl; } }; -- 2.39.5