namespace libstick {
-template<class T>
-std::ostream& operator<<(std::ostream &, const std::vector<T> &);
-
-
/** Persistence computers various persistent homology information on a
* simplex_order of a simplicial_complex. */
template<int MAXDIM, class IT, class VT>
public:
typedef simplicial_complex<MAXDIM, IT, VT> scomplex;
typedef typename scomplex::simplex_order simplex_order;
- typedef typename simplex_order::boundary_matrix boundary_matrix;
+ typedef boolean_colmatrix<IT> boundary_matrix;
typedef boolean_colmatrix<IT> transformation_matrix;
+ typedef std::vector<IT> 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
bm = order.get_boundary_matrix();
rbm = bm;
- tm = create_unit_matrix<transformation_matrix>(bm.size());
-
- reduce_boundary_matrix(rbm, tm);
+ tm = create_unit_matrix<transformation_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;
}
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()) {
os << std::endl;
os << " boundary: " << bm.get_column(c) << std::endl;
- typename boolean_colrowmatrix<IT>::colbase::column_type col = rbm.get_column(c);
+ typename boolean_colmatrix<IT>::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 {
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
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
// 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
}
/** 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;
* invariant. Hence, the resulting b is equal to the product of the boundary
* matrix times v. */
template<class IT, class D>
-void reduce_boundary_matrix(boolean_colrowmatrix<IT> &b, boolean_colmatrix_base<IT, D> &v) {
- assert(b.size() == v.width());
+void reduce_boundary_matrix(boolean_colmatrix_base<IT, D> &b, boolean_colmatrix_base<IT, D> &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<IT> 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<IT>::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<IT>::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<IT>::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<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 << " ";
+ // 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;
}
projects / libstick.git / blobdiff