--- /dev/null
+#ifndef persistence_h_aPeinoXeiFeethuz
+#define persistence_h_aPeinoXeiFeethuz
+
+#include <iostream>
+#include <ostream>
+#include <vector>
+#include <cassert>
+
+#include "simplicialcomplex.h"
+#include "booleanmatrix.h"
+
+
+namespace libstick {
+
+
+/** Persistence computers various persistent homology information on a
+ * simplex_order of a simplicial_complex. */
+template<int MAXDIM, class IT, class VT>
+class persistence {
+
+ public:
+ typedef simplicial_complex<MAXDIM, IT, VT> scomplex;
+ typedef typename scomplex::simplex_order simplex_order;
+ 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
+ * 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 */
+ persistence(const simplex_order &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 simplicial order of this persistence object. */
+ const simplex_order& get_order() const {
+ return order;
+ }
+
+ /** Get boundary matrix 'bm' of 'order', compute reduces boundary
+ * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
+ * tm. */
+ void compute_matrices() {
+ if (done_matrices)
+ return;
+ done_matrices = true;
+
+ bm = order.get_boundary_matrix();
+ rbm = bm;
+ 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).isempty()) {
+ // (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));
+ }
+ }
+
+ // A lowest one remained, recall it
+ if (!rbm.get_column(c).isempty()) {
+ 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 lowestones_matrix& get_lowestones_matrix() const {
+ assert(done_matrices);
+ return lowestones;
+ }
+
+ /** Get the boundary matrix of 'order'. */
+ 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() 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() const {
+ assert(done_matrices);
+ return tm;
+ }
+
+ /** Print the sequence of births and deaths of cycles (homology classes). */
+ void interpret_reduction(std::ostream &os) const {
+ assert(done_matrices);
+
+ for (unsigned c=0; c < bm.width(); ++c) {
+ os << c << ". inserting ";
+
+ switch (bm.get_column(c).size()) {
+ case 0:
+ os << "dummy vertex of dim -1";
+ break;
+ case 1:
+ os << "vertex";
+ break;
+ case 2:
+ os << "edge";
+ break;
+ case 3:
+ os << "triangle";
+ break;
+ case 4:
+ os << "tetrahedron";
+ break;
+ default:
+ os << "simplex";
+ break;
+ }
+
+ os << " \t\tvalue: " << order.get_simplex(c).value << std::endl;
+ os << " boundary: " << bm.get_column(c) << std::endl;
+
+ typename boolean_colmatrix<IT>::column_type col = rbm.get_column(c);
+ if (col.isempty()) {
+ os << " \033[1;32mbirth\033[0;m of a cycle: " << tm.get_column(c) << std::endl;
+ } else {
+ os << " \033[1;31mdeath\033[0;m of a cycle: " << rbm.get_column(c) << std::endl;
+ os << " boundary of: " << tm.get_column(c) << std::endl;
+ }
+ os << std::endl;
+ }
+ }
+
+ /** 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) {
+ // A cycle was born
+ if (rbm.get_column(c).isempty()) {
+ 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[c] > 0)
+ p.death = lowestones[c];
+
+ _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 gives birth or kills, be
+ // or not be, tertium non datur.
+ for (unsigned c=0; c < lowestones.size(); ++c) {
+ // A class died with c
+ if (!rbm.get_column(c).isempty())
+ assert(c == 0 || deaths.count(c) > 0);
+ else
+ assert(births.count(c) > 0);
+ }
+#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) {
+ const IT birth = dia.births[i].birth;
+ const IT death = dia.births[i].death;
+
+ os << " ";
+ if (death > 0) {
+ const scomplex& c = get_order().get_complex();
+ const VT pers = c.simplices[death].value - c.simplices[birth].value;
+ os << birth << "\033[1;31m -> " << death << "\033[0;m" << " \tpers: " << pers;
+ } else {
+ os << "\033[1;32m" << birth << "\033[0;m" << " \t\t[essential]";
+ }
+ os << 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.
+ * 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;
+
+ /** 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];
+};
+
+
+} // namespace libstick
+
+#endif
projects / libstick.git / blobdiff