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