include/libstick/persistence.h

   1 #ifndef persistence_h_aPeinoXeiFeethuz
   2 #define persistence_h_aPeinoXeiFeethuz
   3
   4 #include <iostream>
   5 #include <ostream>
   6 #include <vector>
   7 #include <cassert>
   8
   9 #include "simplicialfunction.h"
  10 #include "booleanmatrix.h"
  11
  12
  13 namespace libstick {
  14
  15
  16 /** Persistence computers various persistent homology information on a
  17  * simplex_order of a simplicial_complex. */
  18 template<int MAXDIM, class IT>
  19 class persistence {
  20
  21     public:
  22         typedef simplicial_complex<MAXDIM, IT> scomplex;
  23         typedef typename scomplex::simplex_order simplex_order;
  24         typedef boolean_colmatrix<IT> boundary_matrix;
  25         typedef boolean_colmatrix<IT> transformation_matrix;
  26         typedef std::vector<IT> lowestones_matrix;
  27
  28         /** A point in a persistence diagram comprises birth- and death-index.
  29          * These are the indices of the simplices of 'order' that give birth
  30          * and death to the class. */
  31         struct diagram_point {
  32             /** 'birth'-th simplex in 'order' gave birth to this cycle. */
  33             IT birth;
  34             /** 'death'-th simplex in 'order' killed this cycle. If 'death' is
  35              * zero, the cycle never dies. Otherwise, 'death' is always larger
  36              * than 'birth'. */
  37             IT death;
  38         };
  39
  40         /** Save the births and deaths of simplices of a fixed dimension, say p. */
  41         struct diagram {
  42             /** Saves the sequence of birth/death-points of p-simplices, sorted
  43              * by birth-index. If death-index is zero, the corresponding
  44              * homology class never dies. */
  45             std::vector<diagram_point> births;
  46
  47
  48             /** Gives the p-th Betti number, i.e., the number of immortal
  49              * p-dimensional homology classes. */
  50             size_t betti() const {
  51                 size_t b = 0;
  52
  53                 for (unsigned i=0; i < births.size(); ++i)
  54                     if (births[i].death == 0)
  55                         b++;
  56
  57                 return b;
  58             }
  59
  60             /** Gives the number of p-dimensional homology classes that are
  61              * born by simplex 'from' or earlier and die after simplex 'to' is
  62              * inserted. */
  63             size_t persistent_betti(IT from, IT to) const {
  64                 assert(from <= to);
  65                 size_t b = 0;
  66
  67                 for (unsigned i=0; i < births.size(); ++i) {
  68                     // All following simplices are born too late.
  69                     if (births[i].birth > from)
  70                         break;
  71                     if (births[i].death > to || births[i].death == 0)
  72                         b++;
  73                 }
  74
  75                 return b;
  76             }
  77         };
  78
  79     public:
  80         /** Create a new peristence object on the given simplex_order */
  81         persistence(const simplex_order &order) :
  82             order(order),
  83             bm(0),
  84             rbm(0),
  85             lowestones(0),
  86             tm(0)
  87             {
  88             reset();
  89         }
  90
  91         /** Reset all results gained from 'order'. */
  92         void reset() {
  93             done_matrices = false;
  94             done_diagrams = false;
  95         }
  96
  97         /** Get simplicial order of this persistence object. */
  98         const simplex_order& get_order() const {
  99             return order;
 100         }
 101
 102         /** Get boundary matrix 'bm' of 'order', compute reduces boundary
 103          * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
 104          * tm. */
 105         void compute_matrices() {
 106             if (done_matrices)
 107                 return;
 108             done_matrices = true;
 109
 110             bm = order.get_boundary_matrix();
 111             rbm = bm;
 112             tm = create_unit_matrix<transformation_matrix>(bm.width());
 113             lowestones = lowestones_matrix(bm.width());
 114
 115             // Make every column reduced, i.e., it is a zero-column or contains only one 1.
 116             for (unsigned c=0; c < rbm.width(); ++c) {
 117                 //if (c % 100 == 0)
 118                 //std::cout << "c = " << c << " (" << (100.0*float(c)/rbm.width()) << " %)" << std::endl;
 119                 // Reduce as long as we need to reduce
 120                 while (!rbm.get_column(c).isempty()) {
 121                     // (r, c) is the lowest one of col
 122                     const IT r = rbm.get_column(c).back();
 123                     // This column contains its lowest one on the same row
 124                     const IT c2 = lowestones[r];
 125                     assert(c2 < c);
 126
 127                     // There is no valid c2, hence we have completely reduced
 128                     if (c2 == 0)
 129                         break;
 130                     // Reduce col by column c2
 131                     else {
 132                         assert(rbm.get(r, c));
 133                         rbm.add_column(c, rbm.get_column(c2));
 134                         tm.add_column(c, tm.get_column(c2));
 135                         assert(!rbm.get(r, c));
 136                     }
 137                 }
 138
 139                 // A lowest one remained, recall it
 140                 if (!rbm.get_column(c).isempty()) {
 141                     const IT r = rbm.get_column(c).back();
 142                     assert(lowestones[r] == 0);
 143                     assert(r < c);
 144                     lowestones[r] = c;
 145                 }
 146             }
 147         }
 148
 149         /** Get the lowest-one matrix of 'order'. */
 150         const lowestones_matrix& get_lowestones_matrix() const {
 151             assert(done_matrices);
 152             return lowestones;
 153         }
 154
 155         /** Get the boundary matrix of 'order'. */
 156         const boundary_matrix& get_boundary_matrix() const {
 157             assert(done_matrices);
 158             return bm;
 159         }
 160
 161         /** Get the reduced boundary matrix of 'order'. */
 162         const boundary_matrix& get_reduced_boundary_matrix() const {
 163             assert(done_matrices);
 164             return rbm;
 165         }
 166
 167         /** Get the transformation matrix of 'order', which transforms the
 168          * boundary matrix to the reduced boundary matrix, when multiplied from
 169          * the right. */
 170         const transformation_matrix& get_transformation_matrix() const {
 171             assert(done_matrices);
 172             return tm;
 173         }
 174
 175         /** Print the sequence of births and deaths of cycles (homology
 176          * classes). If a non-NULL simplicial function is given, also the
 177          * values of the simplices are printed.  */
 178         template<class VT>
 179         void interpret_reduction(std::ostream &os, simplicial_function<MAXDIM, IT, VT> *f) const {
 180             assert(done_matrices);
 181
 182             for (unsigned c=0; c < bm.width(); ++c) {
 183                 os << c << ". inserting ";
 184
 185                 switch (bm.get_column(c).size()) {
 186                     case 0:
 187                         os << "dummy vertex of dim -1";
 188                         break;
 189                     case 1:
 190                         os << "vertex";
 191                         break;
 192                     case 2:
 193                         os << "edge";
 194                         break;
 195                     case 3:
 196                         os << "triangle";
 197                         break;
 198                     case 4:
 199                         os << "tetrahedron";
 200                         break;
 201                     default:
 202                         os << "simplex";
 203                         break;
 204                 }
 205                 os << std::endl;
 206
 207                 if (f != NULL)
 208                     os << "  value: " <<  f->get_value(order.order_to_complex(c)) << std::endl;
 209                 os << "  boundary: " << bm.get_column(c) << std::endl;
 210
 211                 typename boolean_colmatrix<IT>::column_type col = rbm.get_column(c);
 212                 if (col.isempty()) {
 213                     os << "  \033[1;32mbirth\033[0;m of a cycle: " << tm.get_column(c) << std::endl;
 214                 } else {
 215                     os << "  \033[1;31mdeath\033[0;m of a cycle: " << rbm.get_column(c);
 216                     os << ", boundary of: " << tm.get_column(c) << std::endl;
 217                 }
 218                 os << std::endl;
 219             }
 220         }
 221
 222         /** Print the sequence of births and deaths of cycles (homology classes). */
 223         void interpret_reduction(std::ostream &os) const {
 224             interpret_reduction<float>(os, NULL);
 225         }
 226
 227         /** Get the 'dim'-dimensional peristence diagram. */
 228         const diagram& get_persistence_diagram(int dim) const{
 229             assert(-1 <= dim);
 230             assert(dim <= MAXDIM);
 231             assert(done_diagrams);
 232
 233             return diagrams[dim+1];
 234         }
 235
 236         /** Compute persistence diagrams */
 237         void compute_diagrams() {
 238             if (done_diagrams)
 239                 return;
 240             done_diagrams = true;
 241
 242             compute_matrices();
 243
 244 #ifndef NDEBUG
 245             std::set<IT> deaths, births;
 246 #endif
 247             //Clear all diagrams
 248             for (int d=0; d < MAXDIM + 2; ++d)
 249                 diagrams[d] = diagram();
 250
 251             for (unsigned c=0; c < lowestones.size(); ++c) {
 252                 // A cycle was born
 253                 if (rbm.get_column(c).isempty()) {
 254                     const int dim = order.get_simplex(c).dim;
 255
 256                     // Create a diagram point
 257                     diagram_point p;
 258                     p.birth = c;
 259                     p.death = 0;
 260                     // If the class actually dies, get the index of the killing simplex
 261                     if (lowestones[c] > 0)
 262                         p.death = lowestones[c];
 263
 264                     _get_persistence_diagram(dim).births.push_back(p);
 265 #ifndef NDEBUG
 266                     assert(births.count(p.birth) == 0);
 267                     assert(p.death == 0 || deaths.count(p.death) == 0);
 268                     births.insert(p.birth);
 269                     deaths.insert(p.death);
 270 #endif
 271                 }
 272             }
 273 #ifndef NDEBUG
 274             // Shakespeare principle: A simplex either gives birth or kills, be
 275             // or not be, tertium non datur.
 276             for (unsigned c=0; c < lowestones.size(); ++c) {
 277                 // A class died with c
 278                 if (!rbm.get_column(c).isempty())
 279                     assert(c == 0 || deaths.count(c) > 0);
 280                 else
 281                     assert(births.count(c) > 0);
 282             }
 283 #endif
 284         }
 285
 286         /** Print birth and death information in the persistence diagrams. If a
 287          * non-NULL simplicial function is given, also the persistence w.r.t.
 288          * this function is printed. */
 289         template<class VT>
 290         void interpret_persistence(std::ostream &os, simplicial_function<MAXDIM, IT, VT> *f) const {
 291             assert(done_diagrams);
 292
 293             for (int d=-1; d <= MAXDIM; ++d) {
 294                 os << "Dimension " << d << std::endl;
 295                 const diagram &dia = get_persistence_diagram(d);
 296
 297                 for (unsigned i=0; i < dia.births.size(); ++i) {
 298                     const IT birth = dia.births[i].birth;
 299                     const IT death = dia.births[i].death;
 300
 301                     os << "  ";
 302                     if (death > 0) {
 303                         os << birth << "\033[1;31m -> " << death << "\033[0;m";
 304
 305                         if (f != NULL) {
 306                             const VT pers = f->get_value(death) -  f->get_value(birth);
 307                             os << " \tpers: " << pers;
 308                         }
 309                     } else {
 310                         os << "\033[1;32m" << birth << "\033[0;m";
 311                         if (f != NULL)
 312                             os << " \t\t[essential]";
 313                     }
 314                     os << std::endl;
 315                 }
 316             }
 317         }
 318
 319         /** Print birth and death information in the persistence diagrams. */
 320         void interpret_persistence(std::ostream &os) const {
 321             interpret_persistence<float>(os, NULL);
 322         }
 323
 324     private:
 325
 326         diagram& _get_persistence_diagram(int dim) {
 327             assert(-1 <= dim);
 328             assert(dim <= MAXDIM);
 329             assert(done_diagrams);
 330
 331             return diagrams[dim+1];
 332         }
 333
 334         /** The underlying simplex order of this diagram. */
 335         const simplex_order &order;
 336
 337         bool done_matrices;
 338         bool done_diagrams;
 339
 340         /** The boundary matrix of 'order' */
 341         boundary_matrix bm;
 342         /** The reduced boundary matrix of 'order' */
 343         boundary_matrix rbm;
 344         /** Only the lowest ones per column of the reduced boundary matrix.
 345          * lowestones[i] contains the column-index at which the lowest one at
 346          * row i appears, and zero if no such column exists. */
 347         lowestones_matrix lowestones;
 348         /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
 349         transformation_matrix tm;
 350
 351         /** The container for all p-dimensional diagrams. The p-th diagram is
 352          * save at index p+1, with -1 <= p <= MAXDIM. */
 353         diagram diagrams[MAXDIM+2];
 354 };
 355
 356
 357 } // namespace libstick
 358
 359 #endif