From 966ee77918ec9c4c86c2f9d027c86c3510e0e6df Mon Sep 17 00:00:00 2001 From: Stefan Huber Date: Tue, 12 Nov 2013 15:32:09 +0100 Subject: [PATCH] Add computaton of persistence diagrams --- include/libstick-0.1/persistence.h | 200 +++++++++++++++++++++++++++-- tests/persistence.h | 71 ++++++++++ 2 files changed, 257 insertions(+), 14 deletions(-) diff --git a/include/libstick-0.1/persistence.h b/include/libstick-0.1/persistence.h index 937a4cc..c048f13 100644 --- a/include/libstick-0.1/persistence.h +++ b/include/libstick-0.1/persistence.h @@ -4,6 +4,7 @@ #include #include #include +#include #include "simplicialcomplex.h" #include "booleanmatrix.h" @@ -26,14 +27,56 @@ class persistence { typedef typename simplex_order::boundary_matrix boundary_matrix; typedef boolean_colmatrix 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 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 */ @@ -41,14 +84,16 @@ class persistence { 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 @@ -62,35 +107,60 @@ class persistence { bm = order.get_boundary_matrix(); rbm = bm; tm = create_unit_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 "; @@ -129,16 +199,118 @@ class persistence { } } + /** 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 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]; }; diff --git a/tests/persistence.h b/tests/persistence.h index 55774dd..f48c386 100644 --- a/tests/persistence.h +++ b/tests/persistence.h @@ -29,6 +29,7 @@ class persistence_TestSuite: public Test::Suite { otorus(torus) { TEST_ADD(persistence_TestSuite::test_matrix_reduction); + TEST_ADD(persistence_TestSuite::test_betti_numbers); } protected: @@ -164,12 +165,14 @@ class persistence_TestSuite: public Test::Suite { 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(); @@ -193,6 +196,7 @@ class persistence_TestSuite: public Test::Suite { 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(); @@ -211,6 +215,73 @@ class persistence_TestSuite: public Test::Suite { //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); } }; -- 2.39.5