#include <iostream>
#include <ostream>
#include <vector>
+#include <cassert>
#include "simplicialcomplex.h"
#include "booleanmatrix.h"
typedef typename simplex_order::boundary_matrix boundary_matrix;
typedef boolean_colmatrix<IT> transformation_matrix;
- /** A point in a diagram comprises birth- and death-index. These are the
- * indices of the simplices of 'order' that give birth and death to the
- * class. */
+ /** A point in a persistence diagram comprises birth- and death-index.
+ * These are the indices of the simplices of 'order' that give birth
+ * and death to the class. */
struct diagram_point {
+ /** 'birth'-th simplex in 'order' gave birth to this cycle. */
IT birth;
+ /** 'death'-th simplex in 'order' killed this cycle. If 'death' is
+ * zero, the cycle never dies. Otherwise, 'death' is always larger
+ * than 'birth'. */
IT death;
};
+ /** Save the births and deaths of simplices of a fixed dimension, say p. */
+ struct diagram {
+ /** Saves the sequence of birth/death-points of p-simplices, sorted
+ * by birth-index. If death-index is zero, the corresponding
+ * homology class never dies. */
+ std::vector<diagram_point> births;
+
+
+ /** Gives the p-th Betti number, i.e., the number of immortal
+ * p-dimensional homology classes. */
+ size_t betti() const {
+ size_t b = 0;
+
+ for (unsigned i=0; i < births.size(); ++i)
+ if (births[i].death == 0)
+ b++;
+
+ return b;
+ }
+
+ /** Gives the number of p-dimensional homology classes that are
+ * born by simplex 'from' or earlier and die after simplex 'to' is
+ * inserted. */
+ size_t persistent_betti(IT from, IT to) const {
+ assert(from <= to);
+ size_t b = 0;
+
+ for (unsigned i=0; i < births.size(); ++i) {
+ // All following simplices are born too late.
+ if (births[i].birth > from)
+ break;
+ if (births[i].death > to || births[i].death == 0)
+ b++;
+ }
+
+ return b;
+ }
+ };
public:
/** Create a new peristence object on the given simplex_order */
order(order),
bm(0),
rbm(0),
+ lowestones(0),
tm(0)
- {
+ {
reset();
}
/** Reset all results gained from 'order'. */
void reset() {
done_matrices = false;
+ done_diagrams = false;
}
/** Get boundary matrix 'bm' of 'order', compute reduces boundary
bm = order.get_boundary_matrix();
rbm = bm;
tm = create_unit_matrix<transformation_matrix>(bm.size());
- reduce_boundary_matrix(rbm, tm);
+ reduce_boundary_matrix(rbm, tm);
assert(rbm == bm * tm);
+
+ 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);
+
+#ifndef NDEBUG
+ 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) {
+ // 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 IT r = lowestones.get_column(c).back();
+ assert(lowestones.get_column(r).size() == 0);
+ assert(lowestones.get_row(c).size() == 0);
+ }
+ }
+#endif
}
/** Get the boundary matrix of 'order'. */
- const boundary_matrix& get_boundary_matrix() {
- compute_matrices();
+ const boundary_matrix& get_boundary_matrix() const {
+ assert(done_matrices);
return bm;
}
/** Get the reduced boundary matrix of 'order'. */
- const boundary_matrix& get_reduced_boundary_matrix() {
- compute_matrices();
+ const boundary_matrix& get_reduced_boundary_matrix() const {
+ assert(done_matrices);
return rbm;
}
/** Get the transformation matrix of 'order', which transforms the
* boundary matrix to the reduced boundary matrix, when multiplied from
* the right. */
- const transformation_matrix& get_transformation_matrix() {
- compute_matrices();
+ const transformation_matrix& get_transformation_matrix() const {
+ assert(done_matrices);
return tm;
}
-
/** Print the sequence of births and deaths of cycles (homology classes). */
- void interpret_reduction(std::ostream &os) {
- compute_matrices();
+ void interpret_reduction(std::ostream &os) const {
+ assert(done_matrices);
+
for (unsigned c=0; c < bm.size(); ++c) {
os << c << ". inserting ";
}
}
+ /** Get the 'dim'-dimensional peristence diagram. */
+ const diagram& get_persistence_diagram(int dim) const{
+ assert(-1 <= dim);
+ assert(dim <= MAXDIM);
+ assert(done_diagrams);
+
+ return diagrams[dim+1];
+ }
+
+ /** Compute persistence diagrams */
+ void compute_diagrams() {
+ if (done_diagrams)
+ return;
+ done_diagrams = true;
+
+ compute_matrices();
+
+#ifndef NDEBUG
+ std::set<IT> deaths, births;
+#endif
+ //Clear all diagrams
+ for (int d=0; d < MAXDIM + 2; ++d)
+ diagrams[d] = diagram();
+
+ for (unsigned c=0; c < lowestones.size(); ++c) {
+ if (lowestones.get_column(c).size() == 0) {
+ const int dim = order.get_simplex(c).dim;
+
+ // Create a diagram point
+ diagram_point p;
+ 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();
+
+ _get_persistence_diagram(dim).births.push_back(p);
+#ifndef NDEBUG
+ assert(births.count(p.birth) == 0);
+ assert(p.death == 0 || deaths.count(p.death) == 0);
+ births.insert(p.birth);
+ deaths.insert(p.death);
+#endif
+ }
+ }
+#ifndef NDEBUG
+ // Shakespeare principle: A simplex either gaves 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)
+ 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
+ }
+
+ /** Print birth and death information in the persistence diagrams. */
+ void interpret_persistence(std::ostream &os) const {
+ assert(done_diagrams);
+
+ for (int d=-1; d <= MAXDIM; ++d) {
+ os << "Dimension " << d << std::endl;
+ const diagram &dia = get_persistence_diagram(d);
+ for (unsigned i=0; i < dia.births.size(); ++i) {
+ os << " ";
+ if (dia.births[i].death > 0)
+ os << dia.births[i].birth << "\033[1;31m -> " << dia.births[i].death;
+ else
+ os << "\033[1;32m" << dia.births[i].birth;
+ os << "\033[0;m" << std::endl;
+ }
+ }
+ }
+
private:
+ diagram& _get_persistence_diagram(int dim) {
+ assert(-1 <= dim);
+ assert(dim <= MAXDIM);
+ assert(done_diagrams);
+
+ return diagrams[dim+1];
+ }
+
/** The underlying simplex order of this diagram. */
const simplex_order ℴ
bool done_matrices;
+ bool done_diagrams;
+ /** The boundary matrix of 'order' */
boundary_matrix bm;
+ /** 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;
+ /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
transformation_matrix tm;
+
+ /** The container for all p-dimensional diagrams. The p-th diagram is
+ * save at index p+1, with -1 <= p <= MAXDIM. */
+ diagram diagrams[MAXDIM+2];
};
otorus(torus)
{
TEST_ADD(persistence_TestSuite::test_matrix_reduction);
+ TEST_ADD(persistence_TestSuite::test_betti_numbers);
}
protected:
void test_matrix_reduction() {
pers ringp(oring);
+ ringp.compute_matrices();
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);
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();
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();
//std::cout << std::endl;
//ballp.interpret_reduction(std::cout);
+ //torusp.compute_diagrams();
+ //torusp.interpret_persistence(std::cout);
+ //std::cout << std::endl;
+ //std::cout << std::endl;
+ }
+
+ void test_betti_numbers() {
+ pers ringp(oring);
+ ringp.compute_diagrams();
+
+ TEST_ASSERT(ringp.get_persistence_diagram(-1).betti() == 0);
+ TEST_ASSERT(ringp.get_persistence_diagram(0).betti() == 0);
+ TEST_ASSERT(ringp.get_persistence_diagram(1).betti() == 1);
+ TEST_ASSERT(ringp.get_persistence_diagram(2).betti() == 0);
+ TEST_ASSERT(ringp.get_persistence_diagram(3).betti() == 0);
+ for (unsigned d=0; d < 3; ++d)
+ TEST_ASSERT(ringp.get_persistence_diagram(d).persistent_betti(oring.size()-1, oring.size()-1) ==
+ ringp.get_persistence_diagram(d).betti());
+
+
+ pers torusp(otorus);
+ torusp.compute_diagrams();
+
+ TEST_ASSERT(torusp.get_persistence_diagram(-1).betti() == 0);
+ TEST_ASSERT(torusp.get_persistence_diagram(0).betti() == 0);
+ TEST_ASSERT(torusp.get_persistence_diagram(1).betti() == 2);
+ TEST_ASSERT(torusp.get_persistence_diagram(2).betti() == 1);
+ TEST_ASSERT(torusp.get_persistence_diagram(3).betti() == 0);
+ for (unsigned d=0; d < 3; ++d)
+ TEST_ASSERT(torusp.get_persistence_diagram(d).persistent_betti(otorus.size()-1, otorus.size()-1) ==
+ torusp.get_persistence_diagram(d).betti());
+
+ // Torus was made from ring. Hence, persistent betti numbers must correspond to betti numbers of ring.
+ for (unsigned d=0; d < 3; ++d)
+ TEST_ASSERT(torusp.get_persistence_diagram(d).persistent_betti(oring.size()-1, oring.size()-1) ==
+ ringp.get_persistence_diagram(d).betti());
+
+
+ pers ballp(oball);
+ ballp.compute_diagrams();
+
+ //std::cout << std::endl;
+ //ballp.interpret_reduction(std::cout);
+ //std::cout << std::endl;
+
+ TEST_ASSERT(ballp.get_persistence_diagram(-1).betti() == 0);
+ TEST_ASSERT(ballp.get_persistence_diagram(0).betti() == 0);
+ TEST_ASSERT(ballp.get_persistence_diagram(1).betti() == 0);
+ TEST_ASSERT(ballp.get_persistence_diagram(2).betti() == 0);
+ TEST_ASSERT(ballp.get_persistence_diagram(3).betti() == 0);
+ for (unsigned d=0; d < 3; ++d)
+ TEST_ASSERT(ballp.get_persistence_diagram(d).persistent_betti(oball.size()-1, oball.size()-1) ==
+ ballp.get_persistence_diagram(d).betti());
+
+ // Before we inserted the two tetrahedra and the interior triangle, we had a sphere.
+ TEST_ASSERT(ballp.get_persistence_diagram(-1).persistent_betti(oball.size()-4, oball.size()-4) == 0);
+ TEST_ASSERT(ballp.get_persistence_diagram(0).persistent_betti(oball.size()-4, oball.size()-4) == 0);
+ TEST_ASSERT(ballp.get_persistence_diagram(1).persistent_betti(oball.size()-4, oball.size()-4) == 0);
+ TEST_ASSERT(ballp.get_persistence_diagram(2).persistent_betti(oball.size()-4, oball.size()-4) == 1);
+ TEST_ASSERT(ballp.get_persistence_diagram(3).persistent_betti(oball.size()-4, oball.size()-4) == 0);
+
+ // Before we inserted the two tetrahedra, we had a double-sphere.
+ TEST_ASSERT(ballp.get_persistence_diagram(-1).persistent_betti(oball.size()-3, oball.size()-3) == 0);
+ TEST_ASSERT(ballp.get_persistence_diagram(0).persistent_betti(oball.size()-3, oball.size()-3) == 0);
+ TEST_ASSERT(ballp.get_persistence_diagram(1).persistent_betti(oball.size()-3, oball.size()-3) == 0);
+ TEST_ASSERT(ballp.get_persistence_diagram(2).persistent_betti(oball.size()-3, oball.size()-3) == 2);
+ TEST_ASSERT(ballp.get_persistence_diagram(3).persistent_betti(oball.size()-3, oball.size()-3) == 0);
}
};
projects / libstick.git / commitdiff