X-Git-Url: https://git.sthu.org/?p=libstick.git;a=blobdiff_plain;f=include%2Flibstick-0.1%2Fpersistence.h;h=6df70abc7c8345b0f11a86d168fdf3bcd9f4f0f2;hp=0a69be22cb0acf7c974595e4a91a3037aac69c2b;hb=799cd8e53b06b4f9ec7d4dde13ab600717982099;hpb=b3a9951228b854b44c507b32a61e9d36e74709d9 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; }