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