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