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 /** Persistence computers various persistent homology information on a
  17  * simplex_order of a simplicial_complex. */
  18 template<int MAXDIM, class IT, class VT>
  19 class persistence {
  20
  21     public:
  22         typedef simplicial_complex<MAXDIM, IT, VT> 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 boundary matrix 'bm' of 'order', compute reduces boundary
  98          * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
  99          * tm. */
 100         void compute_matrices() {
 101             if (done_matrices)
 102                 return;
 103             done_matrices = true;
 104
 105             bm = order.get_boundary_matrix();
 106             rbm = bm;
 107             tm = create_unit_matrix<transformation_matrix>(bm.width());
 108             lowestones = lowestones_matrix(bm.width());
 109
 110             // Make every column reduced, i.e., it is a zero-column or contains only one 1.
 111             for (unsigned c=0; c < rbm.width(); ++c) {
 112                 //if (c % 100 == 0)
 113                 //std::cout << "c = " << c << " (" << (100.0*float(c)/rbm.width()) << " %)" << std::endl;
 114                 // Reduce as long as we need to reduce
 115                 while (rbm.get_column(c).size() > 0) {
 116                     // (r, c) is the lowest one of col
 117                     const IT r = rbm.get_column(c).back();
 118                     // This column contains its lowest one on the same row
 119                     const IT c2 = lowestones[r];
 120                     assert(c2 < c);
 121
 122                     // There is no valid c2, hence we have completely reduced
 123                     if (c2 == 0)
 124                         break;
 125                     // Reduce col by column c2
 126                     else {
 127                         assert(rbm.get(r, c));
 128                         rbm.add_column(c, rbm.get_column(c2));
 129                         tm.add_column(c, tm.get_column(c2));
 130                         assert(!rbm.get(r, c));
 131                     }
 132                 }
 133
 134                 // A lowest one remained, recall it
 135                 if (rbm.get_column(c).size() > 0) {
 136                     const IT r = rbm.get_column(c).back();
 137                     assert(lowestones[r] == 0);
 138                     assert(r < c);
 139                     lowestones[r] = c;
 140                 }
 141             }
 142         }
 143
 144         /** Get the lowest-one matrix of 'order'. */
 145         const lowestones_matrix& get_lowestones_matrix() const {
 146             assert(done_matrices);
 147             return lowestones;
 148         }
 149
 150         /** Get the boundary matrix of 'order'. */
 151         const boundary_matrix& get_boundary_matrix() const {
 152             assert(done_matrices);
 153             return bm;
 154         }
 155
 156         /** Get the reduced boundary matrix of 'order'. */
 157         const boundary_matrix& get_reduced_boundary_matrix() const {
 158             assert(done_matrices);
 159             return rbm;
 160         }
 161
 162         /** Get the transformation matrix of 'order', which transforms the
 163          * boundary matrix to the reduced boundary matrix, when multiplied from
 164          * the right. */
 165         const transformation_matrix& get_transformation_matrix() const {
 166             assert(done_matrices);
 167             return tm;
 168         }
 169
 170         /** Print the sequence of births and deaths of cycles (homology classes). */
 171         void interpret_reduction(std::ostream &os) const {
 172             assert(done_matrices);
 173
 174             for (unsigned c=0; c < bm.width(); ++c) {
 175                 os << c << ". inserting ";
 176
 177                 switch (bm.get_column(c).size()) {
 178                     case 0:
 179                         os << "dummy vertex of dim -1";
 180                         break;
 181                     case 1:
 182                         os << "vertex";
 183                         break;
 184                     case 2:
 185                         os << "edge";
 186                         break;
 187                     case 3:
 188                         os << "triangle";
 189                         break;
 190                     case 4:
 191                         os << "tetrahedron";
 192                         break;
 193                     default:
 194                         os << "simplex";
 195                         break;
 196                 }
 197
 198                 os << std::endl;
 199                 os << "  boundary: " << bm.get_column(c) << std::endl;
 200
 201                 typename boolean_colmatrix<IT>::column_type col = rbm.get_column(c);
 202                 if (col.size() == 0) {
 203                     os << "  \033[1;32mbirth\033[0;m of a cycle: " << tm.get_column(c) << std::endl;
 204                 } else {
 205                     os << "  \033[1;31mdeath\033[0;m of a cycle: " << rbm.get_column(c) << std::endl;
 206                     os << "  boundary of: " << tm.get_column(c) << std::endl;
 207                 }
 208                 os << std::endl;
 209             }
 210         }
 211
 212         /** Get the 'dim'-dimensional peristence diagram. */
 213         const diagram& get_persistence_diagram(int dim) const{
 214             assert(-1 <= dim);
 215             assert(dim <= MAXDIM);
 216             assert(done_diagrams);
 217
 218             return diagrams[dim+1];
 219         }
 220
 221         /** Compute persistence diagrams */
 222         void compute_diagrams() {
 223             if (done_diagrams)
 224                 return;
 225             done_diagrams = true;
 226
 227             compute_matrices();
 228
 229 #ifndef NDEBUG
 230             std::set<IT> deaths, births;
 231 #endif
 232             //Clear all diagrams
 233             for (int d=0; d < MAXDIM + 2; ++d)
 234                 diagrams[d] = diagram();
 235
 236             for (unsigned c=0; c < lowestones.size(); ++c) {
 237                 // A cycle was born
 238                 if (rbm.get_column(c).size() == 0) {
 239                     const int dim = order.get_simplex(c).dim;
 240
 241                     // Create a diagram point
 242                     diagram_point p;
 243                     p.birth = c;
 244                     p.death = 0;
 245                     // If the class actually dies, get the index of the killing simplex
 246                     if (lowestones[c] > 0)
 247                         p.death = lowestones[c];
 248
 249                     _get_persistence_diagram(dim).births.push_back(p);
 250 #ifndef NDEBUG
 251                     assert(births.count(p.birth) == 0);
 252                     assert(p.death == 0 || deaths.count(p.death) == 0);
 253                     births.insert(p.birth);
 254                     deaths.insert(p.death);
 255 #endif
 256                 }
 257             }
 258 #ifndef NDEBUG
 259             // Shakespeare principle: A simplex either gives birth or kills, be
 260             // or not be, tertium non datur.
 261             for (unsigned c=0; c < lowestones.size(); ++c) {
 262                 // A class died with c
 263                 if (rbm.get_column(c).size() != 0)
 264                     assert(c == 0 || deaths.count(c) > 0);
 265                 else
 266                     assert(births.count(c) > 0);
 267             }
 268 #endif
 269         }
 270
 271         /** Print birth and death information in the persistence diagrams. */
 272         void interpret_persistence(std::ostream &os) const {
 273             assert(done_diagrams);
 274
 275             for (int d=-1; d <= MAXDIM; ++d) {
 276                 os << "Dimension " << d << std::endl;
 277                 const diagram &dia = get_persistence_diagram(d);
 278                 for (unsigned i=0; i < dia.births.size(); ++i) {
 279                     os << "  ";
 280                     if (dia.births[i].death > 0)
 281                         os << dia.births[i].birth << "\033[1;31m -> " << dia.births[i].death;
 282                     else
 283                         os << "\033[1;32m" << dia.births[i].birth;
 284                     os << "\033[0;m" << std::endl;
 285                 }
 286             }
 287         }
 288
 289     private:
 290
 291         diagram& _get_persistence_diagram(int dim) {
 292             assert(-1 <= dim);
 293             assert(dim <= MAXDIM);
 294             assert(done_diagrams);
 295
 296             return diagrams[dim+1];
 297         }
 298
 299         /** The underlying simplex order of this diagram. */
 300         const simplex_order &order;
 301
 302         bool done_matrices;
 303         bool done_diagrams;
 304
 305         /** The boundary matrix of 'order' */
 306         boundary_matrix bm;
 307         /** The reduced boundary matrix of 'order' */
 308         boundary_matrix rbm;
 309         /** Only the lowest ones per column of the reduced boundary matrix.
 310          * lowestones[i] contains the column-index at which the lowest one at
 311          * row i appears, and zero if no such column exists. */
 312         lowestones_matrix lowestones;
 313         /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
 314         transformation_matrix tm;
 315
 316         /** The container for all p-dimensional diagrams. The p-th diagram is
 317          * save at index p+1, with -1 <= p <= MAXDIM. */
 318         diagram diagrams[MAXDIM+2];
 319 };
 320
 321
 322 /** Takes the boundary matrix b and transforms it inplace into a reduced matrix
 323  * b. The matrix v is required to be a unit matrix having the same dimensions
 324  * as b. It is maintained in order to leave the product b*inverse(v) as an
 325  * invariant.  Hence, the resulting b is equal to the product of the boundary
 326  * matrix times v. */
 327 template<class IT, class D>
 328 void reduce_boundary_matrix(boolean_colmatrix_base<IT, D> &b, boolean_colmatrix_base<IT, D> &v) {
 329     assert(b.width() == v.width());
 330
 331     // lowestones[i] gives the column whose lowest 1 is at row i. If it contains
 332     // zero, there is no such column.
 333     std::vector<IT> lowestones(b.width());
 334
 335     // Make every column reduced, i.e., it is a zero-column or contains only one 1.
 336     for (unsigned c=0; c < b.width(); ++c) {
 337         //if (c % 100 == 0)
 338             //std::cout << "c = " << c << " (" << (100.0*float(c)/b.width()) << " %)" << std::endl;
 339         // Reduce as long as we need to reduce
 340         while (b.get_column(c).size() > 0) {
 341             // (r, c) is the lowest one of col
 342             const IT r = b.get_column(c).back();
 343             // This column contains its lowest one on the same row
 344             const IT c2 = lowestones[r];
 345             assert(c2 < c);
 346
 347             // There is no valid c2, hence we have completely reduced
 348             if (c2 == 0)
 349                 break;
 350             // Reduce col by column c2
 351             else {
 352                 assert(b.get(r, c));
 353                 b.add_column(c, b.get_column(c2));
 354                 v.add_column(c, v.get_column(c2));
 355                 assert(!b.get(r, c));
 356             }
 357         }
 358
 359         // A lowest one remained, recall it
 360         if (b.get_column(c).size() > 0) {
 361             const IT r = b.get_column(c).back();
 362             assert(lowestones[r] == 0);
 363             assert(r < c);
 364             lowestones[r] = c;
 365         }
 366     }
 367 }
 368
 369
 370 } // namespace libstick
 371
 372 #endif