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 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 classes). */
 176         void interpret_reduction(std::ostream &os) const {
 177             assert(done_matrices);
 178
 179             for (unsigned c=0; c < bm.width(); ++c) {
 180                 os << c << ". inserting ";
 181
 182                 switch (bm.get_column(c).size()) {
 183                     case 0:
 184                         os << "dummy vertex of dim -1";
 185                         break;
 186                     case 1:
 187                         os << "vertex";
 188                         break;
 189                     case 2:
 190                         os << "edge";
 191                         break;
 192                     case 3:
 193                         os << "triangle";
 194                         break;
 195                     case 4:
 196                         os << "tetrahedron";
 197                         break;
 198                     default:
 199                         os << "simplex";
 200                         break;
 201                 }
 202
 203                 os << " \t\tvalue: " << order.get_simplex(c).value << std::endl;
 204                 os << "  boundary: " << bm.get_column(c) << std::endl;
 205
 206                 typename boolean_colmatrix<IT>::column_type col = rbm.get_column(c);
 207                 if (col.isempty()) {
 208                     os << "  \033[1;32mbirth\033[0;m of a cycle: " << tm.get_column(c) << std::endl;
 209                 } else {
 210                     os << "  \033[1;31mdeath\033[0;m of a cycle: " << rbm.get_column(c) << std::endl;
 211                     os << "  boundary of: " << tm.get_column(c) << std::endl;
 212                 }
 213                 os << std::endl;
 214             }
 215         }
 216
 217         /** Get the 'dim'-dimensional peristence diagram. */
 218         const diagram& get_persistence_diagram(int dim) const{
 219             assert(-1 <= dim);
 220             assert(dim <= MAXDIM);
 221             assert(done_diagrams);
 222
 223             return diagrams[dim+1];
 224         }
 225
 226         /** Compute persistence diagrams */
 227         void compute_diagrams() {
 228             if (done_diagrams)
 229                 return;
 230             done_diagrams = true;
 231
 232             compute_matrices();
 233
 234 #ifndef NDEBUG
 235             std::set<IT> deaths, births;
 236 #endif
 237             //Clear all diagrams
 238             for (int d=0; d < MAXDIM + 2; ++d)
 239                 diagrams[d] = diagram();
 240
 241             for (unsigned c=0; c < lowestones.size(); ++c) {
 242                 // A cycle was born
 243                 if (rbm.get_column(c).isempty()) {
 244                     const int dim = order.get_simplex(c).dim;
 245
 246                     // Create a diagram point
 247                     diagram_point p;
 248                     p.birth = c;
 249                     p.death = 0;
 250                     // If the class actually dies, get the index of the killing simplex
 251                     if (lowestones[c] > 0)
 252                         p.death = lowestones[c];
 253
 254                     _get_persistence_diagram(dim).births.push_back(p);
 255 #ifndef NDEBUG
 256                     assert(births.count(p.birth) == 0);
 257                     assert(p.death == 0 || deaths.count(p.death) == 0);
 258                     births.insert(p.birth);
 259                     deaths.insert(p.death);
 260 #endif
 261                 }
 262             }
 263 #ifndef NDEBUG
 264             // Shakespeare principle: A simplex either gives birth or kills, be
 265             // or not be, tertium non datur.
 266             for (unsigned c=0; c < lowestones.size(); ++c) {
 267                 // A class died with c
 268                 if (!rbm.get_column(c).isempty())
 269                     assert(c == 0 || deaths.count(c) > 0);
 270                 else
 271                     assert(births.count(c) > 0);
 272             }
 273 #endif
 274         }
 275
 276         /** Print birth and death information in the persistence diagrams. */
 277         void interpret_persistence(std::ostream &os) const {
 278             assert(done_diagrams);
 279
 280             for (int d=-1; d <= MAXDIM; ++d) {
 281                 os << "Dimension " << d << std::endl;
 282                 const diagram &dia = get_persistence_diagram(d);
 283
 284                 for (unsigned i=0; i < dia.births.size(); ++i) {
 285                     const IT birth = dia.births[i].birth;
 286                     const IT death = dia.births[i].death;
 287
 288                     os << "  ";
 289                     if (death > 0) {
 290                         const scomplex& c = get_order().get_complex();
 291                         const VT pers = c.simplices[death].value - c.simplices[birth].value;
 292                         os << birth << "\033[1;31m -> " << death << "\033[0;m" << " \tpers: " << pers;
 293                     } else {
 294                         os << "\033[1;32m" << birth << "\033[0;m" << " \t\t[essential]";
 295                     }
 296                     os << std::endl;
 297                 }
 298             }
 299         }
 300
 301     private:
 302
 303         diagram& _get_persistence_diagram(int dim) {
 304             assert(-1 <= dim);
 305             assert(dim <= MAXDIM);
 306             assert(done_diagrams);
 307
 308             return diagrams[dim+1];
 309         }
 310
 311         /** The underlying simplex order of this diagram. */
 312         const simplex_order &order;
 313
 314         bool done_matrices;
 315         bool done_diagrams;
 316
 317         /** The boundary matrix of 'order' */
 318         boundary_matrix bm;
 319         /** The reduced boundary matrix of 'order' */
 320         boundary_matrix rbm;
 321         /** Only the lowest ones per column of the reduced boundary matrix.
 322          * lowestones[i] contains the column-index at which the lowest one at
 323          * row i appears, and zero if no such column exists. */
 324         lowestones_matrix lowestones;
 325         /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
 326         transformation_matrix tm;
 327
 328         /** The container for all p-dimensional diagrams. The p-th diagram is
 329          * save at index p+1, with -1 <= p <= MAXDIM. */
 330         diagram diagrams[MAXDIM+2];
 331 };
 332
 333
 334 } // namespace libstick
 335
 336 #endif