From: Stefan Huber Date: Tue, 19 Nov 2013 15:27:22 +0000 (+0100) Subject: Use more efficient colmatrix-based reduction X-Git-Tag: v0.1~16 X-Git-Url: https://git.sthu.org/?a=commitdiff_plain;h=799cd8e53b06b4f9ec7d4dde13ab600717982099;p=libstick.git Use more efficient colmatrix-based reduction - Use colmatrix as boundary matrix. - Use faster set_symmetric_difference() for colmatrix::add_column(). - Move diagram sanity checks to unit tests. - Make lowestones a col-index vector indexed by row-indices. - Test lowestones in unit tests, which are independt of the reduction algorithm. --- diff --git a/include/libstick-0.1/booleanmatrix.h b/include/libstick-0.1/booleanmatrix.h index 8626f53..2b64a4e 100644 --- a/include/libstick-0.1/booleanmatrix.h +++ b/include/libstick-0.1/booleanmatrix.h @@ -15,6 +15,10 @@ namespace libstick { +template +std::ostream& operator<<(std::ostream &, const std::vector &); + + /** The base class of boolean_colmatrix and boolean_colrowmatrix which implements * the common logic of both. */ template @@ -32,11 +36,26 @@ class boolean_colmatrix_base { } public: - /** Get height resp. width of the matrix. */ + /** A casting constructor for any colmatrix with same entries type. */ + template + boolean_colmatrix_base(const boolean_colmatrix_base &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].size() > 0) + 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()); @@ -61,6 +80,11 @@ class boolean_colmatrix_base { return cols[c]; } + /** Get all columns */ + const std::vector& 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) { @@ -72,12 +96,14 @@ class boolean_colmatrix_base { } /** Two matrices are equal iff they have the same entries */ - bool operator==(const boolean_colmatrix_base &m) const { - return cols == m.cols; + template + bool operator==(const boolean_colmatrix_base &m) const { + return cols == m.get_columns(); } /** Two matrices are equal iff they have the same entries */ - bool operator!=(const boolean_colmatrix_base &m) const { + template + bool operator!=(const boolean_colmatrix_base &m) const { return !(*this == m); } @@ -131,7 +157,7 @@ class boolean_colmatrix_base { #endif } - private: + protected: /** The matrix is the set of columns. */ std::vector cols; @@ -151,6 +177,7 @@ class boolean_colmatrix : public boolean_colmatrix_base > base; + typedef typename base::column_type column_type; /** Create a matrix with 'width' columns, initalized with zero entries. */ boolean_colmatrix(size_t columns) : @@ -161,6 +188,37 @@ class boolean_colmatrix : public boolean_colmatrix_base(this); } @@ -314,6 +372,23 @@ class boolean_colrowmatrix : public boolean_colmatrix_base + boolean_colrowmatrix(const boolean_colmatrix_base& 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[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[i]) == true); + } + } + /** Override implementation. */ void _set(index_type r, index_type c, bool value) { colbase::_set(r, c, value); @@ -352,6 +427,12 @@ class boolean_colrowmatrix : public boolean_colmatrix_base + bool operator==(const boolean_colmatrix_base &m) const { + return colbase::operator==(m); + } + /** Two matrices are equal iff they have the same entries */ bool operator!=(const boolean_colrowmatrix &m) const { return !(*this == m); @@ -524,6 +605,21 @@ std::ostream& operator<<(std::ostream &os, boolean_colmatrix_base &mat) { return os << m; } +template +std::ostream& operator<<(std::ostream& os, const std::vector &vec) { + os << "["; + + typename std::vector::const_iterator it = vec.begin(); + while ( it != vec.end()) { + os << *it; + if (++it != vec.end()) + os << " "; + } + + os << "]"; + return os; +} + } #endif diff --git a/include/libstick-0.1/persistence.h b/include/libstick-0.1/persistence.h index 0a69be2..6df70ab 100644 --- a/include/libstick-0.1/persistence.h +++ b/include/libstick-0.1/persistence.h @@ -13,10 +13,6 @@ namespace libstick { -template -std::ostream& operator<<(std::ostream &, const std::vector &); - - /** Persistence computers various persistent homology information on a * simplex_order of a simplicial_complex. */ template @@ -25,8 +21,9 @@ class persistence { public: typedef simplicial_complex scomplex; typedef typename scomplex::simplex_order simplex_order; - typedef typename simplex_order::boundary_matrix boundary_matrix; + typedef boolean_colmatrix boundary_matrix; typedef boolean_colmatrix transformation_matrix; + typedef std::vector lowestones_matrix; /** A point in a persistence diagram comprises birth- and death-index. * These are the indices of the simplices of 'order' that give birth @@ -107,18 +104,45 @@ class persistence { bm = order.get_boundary_matrix(); rbm = bm; - tm = create_unit_matrix(bm.size()); - - reduce_boundary_matrix(rbm, tm); + tm = create_unit_matrix(bm.width()); + lowestones = lowestones_matrix(bm.width()); + + // Make every column reduced, i.e., it is a zero-column or contains only one 1. + for (unsigned c=0; c < rbm.width(); ++c) { + //if (c % 100 == 0) + //std::cout << "c = " << c << " (" << (100.0*float(c)/rbm.width()) << " %)" << std::endl; + // Reduce as long as we need to reduce + while (rbm.get_column(c).size() > 0) { + // (r, c) is the lowest one of col + const IT r = rbm.get_column(c).back(); + // This column contains its lowest one on the same row + const IT c2 = lowestones[r]; + assert(c2 < c); + + // There is no valid c2, hence we have completely reduced + if (c2 == 0) + break; + // Reduce col by column c2 + else { + assert(rbm.get(r, c)); + rbm.add_column(c, rbm.get_column(c2)); + tm.add_column(c, tm.get_column(c2)); + assert(!rbm.get(r, c)); + } + } - lowestones = boundary_matrix(rbm.size()); - for (unsigned c=0; c < rbm.size(); ++c) - if (rbm.get_column(c).size() > 0) - lowestones.set(rbm.get_column(c).back(), c, true); + // A lowest one remained, recall it + if (rbm.get_column(c).size() > 0) { + const IT r = rbm.get_column(c).back(); + assert(lowestones[r] == 0); + assert(r < c); + lowestones[r] = c; + } + } } /** Get the lowest-one matrix of 'order'. */ - const boundary_matrix& get_lowestones_matrix() const { + const lowestones_matrix& get_lowestones_matrix() const { assert(done_matrices); return lowestones; } @@ -147,7 +171,7 @@ class persistence { void interpret_reduction(std::ostream &os) const { assert(done_matrices); - for (unsigned c=0; c < bm.size(); ++c) { + for (unsigned c=0; c < bm.width(); ++c) { os << c << ". inserting "; switch (bm.get_column(c).size()) { @@ -174,7 +198,7 @@ class persistence { os << std::endl; os << " boundary: " << bm.get_column(c) << std::endl; - typename boolean_colrowmatrix::colbase::column_type col = rbm.get_column(c); + typename boolean_colmatrix::column_type col = rbm.get_column(c); if (col.size() == 0) { os << " \033[1;32mbirth\033[0;m of a cycle: " << tm.get_column(c) << std::endl; } else { @@ -210,7 +234,8 @@ class persistence { diagrams[d] = diagram(); for (unsigned c=0; c < lowestones.size(); ++c) { - if (lowestones.get_column(c).size() == 0) { + // A cycle was born + if (rbm.get_column(c).size() == 0) { const int dim = order.get_simplex(c).dim; // Create a diagram point @@ -218,8 +243,8 @@ class persistence { p.birth = c; p.death = 0; // If the class actually dies, get the index of the killing simplex - if (lowestones.get_row(c).size() > 0) - p.death = lowestones.get_row(c).back(); + if (lowestones[c] > 0) + p.death = lowestones[c]; _get_persistence_diagram(dim).births.push_back(p); #ifndef NDEBUG @@ -234,18 +259,12 @@ class persistence { // Shakespeare principle: A simplex either gives birth or kills, be // or not be, tertium non datur. for (unsigned c=0; c < lowestones.size(); ++c) { - if (lowestones.get_column(c).size() != 0) + // A class died with c + if (rbm.get_column(c).size() != 0) assert(c == 0 || deaths.count(c) > 0); else assert(births.count(c) > 0); } - // Even Jesus died after his birth. Homology classes are nothing better. - for (int d=-1; d <= MAXDIM; ++d) { - const diagram &dia = get_persistence_diagram(d); - for (unsigned i=0; i < dia.births.size(); ++i) - if (dia.births[i].death > 0) - assert(dia.births[i].birth < dia.births[i].death); - } #endif } @@ -288,9 +307,9 @@ class persistence { /** The reduced boundary matrix of 'order' */ boundary_matrix rbm; /** Only the lowest ones per column of the reduced boundary matrix. - * Each one at (r, c) corresponds to the death of a class when - * inserting simplex c, which was born when inserting simplex. */ - boundary_matrix lowestones; + * lowestones[i] contains the column-index at which the lowest one at + * row i appears, and zero if no such column exists. */ + lowestones_matrix lowestones; /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */ transformation_matrix tm; @@ -306,51 +325,45 @@ class persistence { * invariant. Hence, the resulting b is equal to the product of the boundary * matrix times v. */ template -void reduce_boundary_matrix(boolean_colrowmatrix &b, boolean_colmatrix_base &v) { - assert(b.size() == v.width()); +void reduce_boundary_matrix(boolean_colmatrix_base &b, boolean_colmatrix_base &v) { + assert(b.width() == v.width()); + + // lowestones[i] gives the column whose lowest 1 is at row i. If it contains + // zero, there is no such column. + std::vector lowestones(b.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 boolean_colrowmatrix::colbase::column_type &col = b.get_column(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 boolean_colrowmatrix::rowbase::row_type row(b.get_row(r).size()); - copy(b.get_row(r).begin(), b.get_row(r).end(), row.begin()); - assert(row.size() >= 1); - - // Get rid of 1s at that row right of column c. - typename boolean_colrowmatrix::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.add_column(*it, col); - v.add_column(*it, v.get_column(c)); - assert(!b.get(r, *it)); + for (unsigned c=0; c < b.width(); ++c) { + //if (c % 100 == 0) + //std::cout << "c = " << c << " (" << (100.0*float(c)/b.width()) << " %)" << std::endl; + // Reduce as long as we need to reduce + while (b.get_column(c).size() > 0) { + // (r, c) is the lowest one of col + const IT r = b.get_column(c).back(); + // This column contains its lowest one on the same row + const IT c2 = lowestones[r]; + assert(c2 < c); + + // There is no valid c2, hence we have completely reduced + if (c2 == 0) + break; + // Reduce col by column c2 + else { + assert(b.get(r, c)); + b.add_column(c, b.get_column(c2)); + v.add_column(c, v.get_column(c2)); + assert(!b.get(r, c)); + } } - } -} - -template -std::ostream& operator<<(std::ostream& os, const std::vector &vec) { - os << "["; - typename std::vector::const_iterator it = vec.begin(); - while ( it != vec.end()) { - os << *it; - if (++it != vec.end()) - os << " "; + // A lowest one remained, recall it + if (b.get_column(c).size() > 0) { + const IT r = b.get_column(c).back(); + assert(lowestones[r] == 0); + assert(r < c); + lowestones[r] = c; + } } - - os << "]"; - return os; } diff --git a/include/libstick-0.1/simplicialcomplex.h b/include/libstick-0.1/simplicialcomplex.h index 1183a92..973d632 100644 --- a/include/libstick-0.1/simplicialcomplex.h +++ b/include/libstick-0.1/simplicialcomplex.h @@ -92,7 +92,7 @@ class simplicial_complex { class simplex_order { public: - typedef boolean_colrowmatrix boundary_matrix; + typedef boolean_colmatrix boundary_matrix; /** Create a standard order of the complex c, i.e., the identity permutation. */ simplex_order(const simplcompltype &c) : diff --git a/tests/booleanmatrix.h b/tests/booleanmatrix.h index 71533ff..6cbf538 100644 --- a/tests/booleanmatrix.h +++ b/tests/booleanmatrix.h @@ -63,10 +63,10 @@ class boolean_matrix_TestSuite: public Test::Suite { for (unsigned c=0; c < size; ++c) TEST_ASSERT(mat.get_column(c).size() == 0); + col.push_back(0); col.push_back(3); - col.push_back(14); col.push_back(8); - col.push_back(0); + col.push_back(14); // Add column and test for values mat.add_column(5, col); diff --git a/tests/image.h b/tests/image.h index ff7817a..8633a56 100644 --- a/tests/image.h +++ b/tests/image.h @@ -49,7 +49,7 @@ class image_TestSuite: public Test::Suite { sorder so(s); assert(so.is_filtration()); - typedef sorder::boundary_matrix bmatrix; + typedef boolean_colrowmatrix bmatrix; bmatrix bm = so.get_boundary_matrix(); // Check for the right vertex incidences diff --git a/tests/persistence.h b/tests/persistence.h index 27b01f2..832457e 100644 --- a/tests/persistence.h +++ b/tests/persistence.h @@ -15,6 +15,7 @@ class persistence_TestSuite: public Test::Suite { typedef simplicial_complex<3, uint32_t, double> scomplex; typedef persistence<3, uint32_t, double> pers; typedef pers::boundary_matrix bm; + typedef pers::lowestones_matrix lom; typedef pers::transformation_matrix tm; bool setupcalled; @@ -31,6 +32,7 @@ class persistence_TestSuite: public Test::Suite { TEST_ADD(persistence_TestSuite::test_matrix_reduction); TEST_ADD(persistence_TestSuite::test_betti_numbers); TEST_ADD(persistence_TestSuite::test_lowestones); + TEST_ADD(persistence_TestSuite::test_diagram); } protected: @@ -170,48 +172,37 @@ class persistence_TestSuite: public Test::Suite { const bm &ringbm = ringp.get_boundary_matrix(); const bm &ringrbm = ringp.get_reduced_boundary_matrix(); const tm &ringtm = ringp.get_transformation_matrix(); - TEST_ASSERT(ringrbm == ringbm * ringtm); + TEST_ASSERT(ringrbm == boolean_colrowmatrix(ringbm) * ringtm); + pers torusp(otorus); torusp.compute_matrices(); const bm &torusbm = torusp.get_boundary_matrix(); const bm &torusrbm = torusp.get_reduced_boundary_matrix(); const tm &torustm = torusp.get_transformation_matrix(); - TEST_ASSERT(torusrbm == torusbm * torustm); - - uint32_t torusrbme_coords[][2] = { - {0, 1}, {2, 7}, {3, 7}, {1, 8}, {2, 8}, {4, 10}, {5, 10}, {4, 11}, {6, 11}, {1, 13}, - {4, 13}, {10, 15}, {13, 15}, {14, 15}, {9, 19}, {10, 19}, {13, 19}, {17, 19}, {7, 20}, - {16, 20}, {18, 20}, {7, 21}, {9, 21}, {10, 21}, {11, 21}, {13, 21}, {16, 21}, {12, 23}, - {13, 23}, {22, 23}, {7, 24}, {8, 24}, {9, 24}, {10, 24}, {11, 24}, {12, 24}, {25, 28}, - {26, 28}, {25, 29}, {27, 29}, {1, 31}, {25, 31}, {28, 33}, {31, 33}, {32, 33}, {30, 35}, - {31, 35}, {34, 35}, {8, 37}, {30, 37}, {31, 37}, {36, 37}, {9, 40}, {28, 40}, {31, 40}, - {39, 40}, {9, 41}, {28, 41}, {29, 41}, {31, 41}, {38, 41}, {7, 42}, {8, 42}, {9, 42}, - {28, 42}, {29, 42}, {30, 42}, {7, 49}, {8, 49}, {9, 49}, {28, 49}, {29, 49}, {43, 49}, - {47, 49}, {10, 50}, {11, 50}, {28, 50}, {29, 50}, {43, 50}, {45, 50}, {29, 51}, {44, 51}, - {48, 51}, {10, 52}, {28, 52}, {43, 52}, {44, 52}, {28, 53}, {43, 53}, {46, 53 } - }; - bm torusrbme(torusbm.size()); - torusrbme.set_all(torusrbme_coords, sizeof(torusrbme_coords)/(2*sizeof(uint32_t)), true); - TEST_ASSERT(torusrbme == torusrbm); + TEST_ASSERT(torusrbm == boolean_colrowmatrix(torusbm) * torustm); + + const lom &toruslo = torusp.get_lowestones_matrix(); + uint32_t torusloe_coords[] = {1, 0, 8, 7, 13, 10, 11, 0, 0, 0, 0, + 0, 24, 0, 15, 0, 21, 19, 20, 0, 0, 0, 23, 0, 0, 31, 28, 29, 0, + 0, 42, 0, 33, 0, 35, 0, 37, 0, 41, 40, 0, 0, 0, 0, 52, 50, 53, + 49, 51, 0, 0, 0, 0, 0, 0}; + lom torusloe(torusloe_coords, torusloe_coords + sizeof(torusloe_coords)/sizeof(uint32_t)); + TEST_ASSERT(torusloe == toruslo); + pers ballp(oball); ballp.compute_matrices(); const bm &ballbm = ballp.get_boundary_matrix(); const bm &ballrbm = ballp.get_reduced_boundary_matrix(); const tm &balltm = ballp.get_transformation_matrix(); - TEST_ASSERT(ballrbm == ballbm * balltm); - - uint32_t ballrbme_coords[][2] = { - {0, 1}, {1, 6}, {2, 6}, {1, 7}, {3, 7}, {1, 8}, {4, 8}, {1, 10}, {5, 10}, {6, 14}, - {10, 14}, {11, 14}, {6, 15}, {7, 15}, {10, 15}, {12, 15}, {6, 16}, {7, 16}, {8, 16}, - {10, 16}, {13, 16}, {6, 17}, {7, 17}, {8, 17}, {9, 17}, {6, 19}, {7, 19}, {18, 19}, - {14, 22}, {15, 22}, {19, 22}, {21, 22}, {14, 23}, {15, 23}, {16, 23}, {17, 23}, {19, 23}, - {20, 23 } - }; - bm ballrbme(ballrbm.size()); - ballrbme.set_all(ballrbme_coords, sizeof(ballrbme_coords)/(2*sizeof(uint32_t)), true); - TEST_ASSERT(ballrbme == ballrbm); + TEST_ASSERT(ballrbm == boolean_colrowmatrix(ballbm) * balltm); + + const lom &balllo = ballp.get_lowestones_matrix(); + uint32_t ballloe_coords[] = {1, 0, 6, 7, 8, 10, 0, 0, 0, 17, 0, 14, + 15, 16, 0, 0, 0, 0, 19, 0, 23, 22, 0, 0}; + lom ballloe(ballloe_coords, ballloe_coords + sizeof(ballloe_coords)/sizeof(uint32_t)); + TEST_ASSERT(ballloe == balllo); //std::cout << std::endl; //ballp.interpret_reduction(std::cout); @@ -293,26 +284,56 @@ class persistence_TestSuite: public Test::Suite { void test_lowestones_impl(const scomplex::simplex_order& so) { pers p(so); p.compute_matrices(); - const bm &lowestones = p.get_lowestones_matrix(); - - for (unsigned c=0; c < lowestones.size(); ++c) - assert(lowestones.get_column(c).size() <= 1); - for (unsigned r=0; r < lowestones.size(); ++r) - assert(lowestones.get_row(r).size() <= 1); - for (unsigned c=0; c < lowestones.size(); ++c) { + const lom &lowestones = p.get_lowestones_matrix(); + + // Check that lowestones contains every column-index only once, + // except the zeros. + lom cpy = lowestones; + sort(cpy.begin(), cpy.end()); + lom::iterator it = cpy.begin(); + // Skip the zeros + while (it != cpy.end() && *it == 0) + ++it; + // From now on, no column index appears twice + for (; it+1 != cpy.end(); ++it) + TEST_ASSERT(*it < *(it+1)); + + for (unsigned r=0; r < lowestones.size(); ++r) { // If (r,c) is a one then // (i) a cycle dies with c --> row c is zero // (ii) a cycle is born with r --> column r is zero //Hence // (i) the column r is a zero-column // (i) the row c is a zero-column - if (lowestones.get_column(c).size() > 0) { - const uint32_t r = lowestones.get_column(c).back(); - assert(lowestones.get_column(r).size() == 0); - assert(lowestones.get_row(c).size() == 0); + + const uint32_t c = lowestones[r]; + if (c > 0) { + TEST_ASSERT(p.get_reduced_boundary_matrix().get(r,c) == true); + TEST_ASSERT(p.get_reduced_boundary_matrix().get_column(r).size() == 0); + TEST_ASSERT(lowestones[c] == 0); } } } + + void test_diagram() { + test_diagram_impl(oring); + test_diagram_impl(otorus); + test_diagram_impl(oball); + } + + void test_diagram_impl(const scomplex::simplex_order& so) { + pers p(so); + p.compute_diagrams(); + + // Even Jesus died after his birth. Homology classes are nothing better. + for (int d=-1; d <= 3; ++d) { + const pers::diagram &dia = p.get_persistence_diagram(d); + for (unsigned i=0; i < dia.births.size(); ++i) + if (dia.births[i].death > 0) + assert(dia.births[i].birth < dia.births[i].death); + } + + } }; #endif