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
   8 #include "simplicialcomplex.h"
   9 #include "booleanmatrix.h"
  10
  11
  12 namespace libstick {
  13
  14
  15 template<class T>
  16 std::ostream& operator<<(std::ostream &, const std::vector<T> &);
  17
  18
  19 /** Persistence computers various persistent homology information on a
  20  * simplex_order of a simplicial_complex. */
  21 template<int MAXDIM, class IT, class VT>
  22 class persistence {
  23
  24     public:
  25         typedef typename simplicial_complex<MAXDIM, IT, VT>::simplex_order simplex_order;
  26         typedef typename simplex_order::boundary_matrix boundary_matrix;
  27         typedef boolean_colmatrix<IT> transformation_matrix;
  28
  29         /** A point in a diagram comprises birth- and death-index. These are the
  30          * indices of the simplices of 'order' that give birth and death to the
  31          * class. */
  32         struct diagram_point {
  33             IT birth;
  34             IT death;
  35         };
  36
  37
  38     public:
  39         /** Create a new peristence object on the given simplex_order */
  40         persistence(const simplex_order &order) :
  41             order(order),
  42             bm(0),
  43             rbm(0),
  44             tm(0)
  45         {
  46             reset();
  47         }
  48
  49         /** Reset all results gained from 'order'. */
  50         void reset() {
  51             done_matrices = false;
  52         }
  53
  54         /** Get boundary matrix 'bm' of 'order', compute reduces boundary
  55          * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
  56          * tm. */
  57         void compute_matrices() {
  58             if (done_matrices)
  59                 return;
  60             done_matrices = true;
  61
  62             bm = order.get_boundary_matrix();
  63             rbm = bm;
  64             tm = create_unit_matrix<transformation_matrix>(bm.size());
  65             reduce_boundary_matrix(rbm, tm);
  66
  67             assert(rbm == bm * tm);
  68         }
  69
  70         /** Get the boundary matrix of 'order'. */
  71         const boundary_matrix& get_boundary_matrix() {
  72             compute_matrices();
  73             return bm;
  74         }
  75
  76         /** Get the reduced boundary matrix of 'order'. */
  77         const boundary_matrix& get_reduced_boundary_matrix() {
  78             compute_matrices();
  79             return rbm;
  80         }
  81
  82         /** Get the transformation matrix of 'order', which transforms the
  83          * boundary matrix to the reduced boundary matrix, when multiplied from
  84          * the right. */
  85         const transformation_matrix& get_transformation_matrix() {
  86             compute_matrices();
  87             return tm;
  88         }
  89
  90
  91         /** Print the sequence of births and deaths of cycles (homology classes). */
  92         void interpret_reduction(std::ostream &os) {
  93             compute_matrices();
  94             for (unsigned c=0; c < bm.size(); ++c) {
  95                 os << c << ". inserting ";
  96
  97                 switch (bm.get_column(c).size()) {
  98                     case 0:
  99                         os << "dummy vertex of dim -1";
 100                         break;
 101                     case 1:
 102                         os << "vertex";
 103                         break;
 104                     case 2:
 105                         os << "edge";
 106                         break;
 107                     case 3:
 108                         os << "triangle";
 109                         break;
 110                     case 4:
 111                         os << "tetrahedron";
 112                         break;
 113                     default:
 114                         os << "simplex";
 115                         break;
 116                 }
 117
 118                 os << std::endl;
 119                 os << "  boundary: " << bm.get_column(c) << std::endl;
 120
 121                 typename boolean_colrowmatrix<IT>::colbase::column_type col = rbm.get_column(c);
 122                 if (col.size() == 0) {
 123                     os << "  \033[1;32mbirth\033[0;m of a cycle: " << tm.get_column(c) << std::endl;
 124                 } else {
 125                     os << "  \033[1;31mdeath\033[0;m of a cycle: " << rbm.get_column(c) << std::endl;
 126                     os << "  boundary of: " << tm.get_column(c) << std::endl;
 127                 }
 128                 os << std::endl;
 129             }
 130         }
 131
 132     private:
 133
 134         /** The underlying simplex order of this diagram. */
 135         const simplex_order &order;
 136
 137         bool done_matrices;
 138
 139         boundary_matrix bm;
 140         boundary_matrix rbm;
 141         transformation_matrix tm;
 142 };
 143
 144
 145 /** Takes the boundary matrix b and transforms it inplace into a reduced matrix
 146  * b. The matrix v is required to be a unit matrix having the same dimensions
 147  * as b. It is maintained in order to leave the product b*inverse(v) as an
 148  * invariant.  Hence, the resulting b is equal to the product of the boundary
 149  * matrix times v. */
 150 template<class IT, class D>
 151 void reduce_boundary_matrix(boolean_colrowmatrix<IT> &b, boolean_colmatrix_base<IT, D> &v) {
 152     assert(b.size() == v.width());
 153
 154     // Make every column reduced, i.e., it is a zero-column or contains only one 1.
 155     for (unsigned c=0; c < b.size(); ++c) {
 156         const typename boolean_colrowmatrix<IT>::colbase::column_type &col = b.get_column(c);
 157         if (col.size() == 0)
 158             continue;
 159
 160         // The column is not a zero-column, take the lowest 1, and reduce the other columns at the row
 161         IT r = col.back();
 162         assert(b.get(r, c));
 163
 164         // Get a copy of the r-th row
 165         typename boolean_colrowmatrix<IT>::rowbase::row_type row(b.get_row(r).size());
 166         copy(b.get_row(r).begin(), b.get_row(r).end(), row.begin());
 167         assert(row.size() >= 1);
 168
 169         // Get rid of 1s at that row right of column c.
 170         typename boolean_colrowmatrix<IT>::row_type::const_iterator it = lower_bound(row.begin(), row.end(), c);
 171         for(; it != row.end(); ++it) {
 172             if (*it <= c)
 173                 continue;
 174
 175             assert(b.get(r, *it));
 176             b.add_column(*it, col);
 177             v.add_column(*it, v.get_column(c));
 178             assert(!b.get(r, *it));
 179         }
 180     }
 181 }
 182
 183 template<class T>
 184 std::ostream& operator<<(std::ostream& os, const std::vector<T> &vec) {
 185     os << "[";
 186
 187     typename std::vector<T>::const_iterator it = vec.begin();
 188     while ( it != vec.end()) {
 189         os << *it;
 190         if (++it != vec.end())
 191             os << " ";
 192     }
 193
 194     os << "]";
 195     return os;
 196 }
 197
 198
 199 } // namespace libstick
 200
 201 #endif