+ /** 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;
+ }
+ }
+ }
+