include/libstick/simplicialcomplex.h

   1 #ifndef simplicialcomplex_h_nealaezeojeeChuh
   2 #define simplicialcomplex_h_nealaezeojeeChuh
   3
   4 #include <stdint.h>
   5 #include <cstdlib>
   6 #include <cstring>
   7 #include <algorithm>
   8 #include <vector>
   9 #include <limits>
  10
  11 #include <iostream>
  12
  13 #include "booleanmatrix.h"
  14
  15
  16 namespace libstick {
  17
  18 /** A simplicial complex is a collection of simplices such that each face is
  19  * also part of the complex. Every simplex has dimension at most MAXDIM. The
  20  * indices of simplices and their faces are of type IT. When a
  21  * simplicial_complex is instantiated, a single (-1) dimensional simplex is
  22  * automatically created. (When we compute persistent homology, we actually
  23  * compute reduced homology.) Each 0-dimensional simplex automatically has this
  24  * simplex as its face.
  25  *
  26  * Based on a simplicial complex, we define a simplicial order. An order
  27  * can be seen as a permutation of the simplices of the complex. If the order
  28  * has the property that every prefix of the permutation is again a complex,
  29  * i.e., all faces of all simplices of each prefix are contained in the prefix,
  30  * then we call it a filtration. The innner class simplex_order gives the
  31  * extended boundary matrix, which can then be used in the persistence class to
  32  * compute persistent homology.
  33  *
  34  * One way to obtain a filtration is to define a monotone simplicial function
  35  * on the complex. That is, each simplex gets a value of type VT assigned. Then
  36  * one can obtain a sub-level set filtration from the complex w.r.t. to this
  37  * simplicial function.
  38  */
  39 template<int MAXDIM, class IT>
  40 class simplicial_complex {
  41
  42     public:
  43         /** The type of this class. */
  44         typedef simplicial_complex<MAXDIM, IT> scomplex;
  45         /** Type of indices of simplices. */
  46         typedef IT index_type;
  47
  48         /** A simplex of the complex. This struct is an aggregate. */
  49         struct simplex {
  50             /** Dimension of the simplex. */
  51             int dim;
  52             /** The indices of the faces of the simplex. */
  53             index_type faces[MAXDIM+1];
  54
  55             /** Create a new simplex with dimension 'dim' and (dim+1)-faces. If
  56              * the simplex is 0-dimensional, its face is automatically set to
  57              * one (-1)-dimensional simplex. */
  58             static simplex create_simplex(int dim, index_type* faces) {
  59                 assert(0 <= dim && dim <= MAXDIM);
  60
  61                 simplex s = {dim, {0}};
  62
  63                 if (dim >= 1)
  64                     memcpy(s.faces, faces, face_count_bydim(dim)*sizeof(index_type));
  65
  66                 return s;
  67             }
  68
  69             /** Create a (-1)-dimensional simplex. */
  70             static simplex create_minusonedim_simplex() {
  71                 simplex s = {-1, {0}};
  72                 return s;
  73             }
  74
  75             /** Get number of faces. */
  76             size_t face_count() const {
  77                 return face_count_bydim(dim);
  78             }
  79
  80             /** Get number of faces of a dim-dimensional simplex. */
  81             static size_t face_count_bydim(int dim) {
  82                 assert(-1 <= dim && dim <= MAXDIM);
  83                 return dim + 1;
  84             }
  85         };
  86
  87         /** An order of the simplices of complex c. An order can be interpreted
  88          * as a permuation of the complex's std::vector of simplices.  */
  89         class simplex_order {
  90
  91             public:
  92                 typedef boolean_colmatrix<IT> boundary_matrix;
  93                 typedef std::vector<IT> order_type;
  94
  95                 /** Create a standard order of the complex c, i.e., the identity permutation. */
  96                 simplex_order(const scomplex &c) :
  97                     c(c) {
  98                     reset();
  99                 }
 100
 101                 /** Reset order to the identity permutation of the complex's simplices. */
 102                 void reset() {
 103                     order.clear();
 104                     for (unsigned i=0; i < c.size(); ++i)
 105                         order.push_back(i);
 106                     revorder = order;
 107                 }
 108
 109                 /** Return number of simplices. */
 110                 size_t size() const {
 111                     assert(order.size() == revorder.size());
 112                     return order.size();
 113                 }
 114
 115                 /** Get i-th simplex in the simplex order. */
 116                 const simplex& get_simplex(index_type i) const {
 117                     return c.simplices.at(order_to_complex(i));
 118                 }
 119
 120                 /** Translate an index in the order to a simplex in the complex. */
 121                 index_type order_to_complex(index_type i) const {
 122                     assert(0 <= i && i < size());
 123                     return order.at(i);
 124                 }
 125
 126                 /** Translate an index in the complex to a simplex in the order. */
 127                 index_type complex_to_order(index_type i) const {
 128                     assert(0 <= i && i < size());
 129                     return revorder.at(i);
 130                 }
 131
 132                 /** Return the underlying complex. */
 133                 const scomplex& get_complex() const {
 134                     return c;
 135                 }
 136
 137                 /** Check if this order is a permutation of the underlying
 138                  * complex. */
 139                 bool is_permutation() const {
 140                     if (size() != c.size())
 141                         return false;
 142
 143                     // Sort the index array and check if every index is present.
 144                     order_type o = order;
 145                     std::sort(o.begin(), o.end());
 146                     for (unsigned i=0; i < size(); ++i)
 147                         if (o.at(i) != i)
 148                             return false;
 149
 150                     return true;
 151                 }
 152
 153                 /** Returns true iff the faces of simplex i are before i in
 154                  * this order. */
 155                 bool is_filtration() const {
 156                     assert(size() == c.size());
 157
 158                     for (unsigned i=0; i < size(); ++i)
 159                         for (unsigned f=0; f < get_simplex(i).face_count(); ++f)
 160                             if (complex_to_order(get_simplex(i).faces[f]) >= i)
 161                                 return false;
 162
 163                     return true;
 164                 }
 165
 166                 /** Randomize order. It has hardly any impact on runtime, but
 167                  * it makes cycles "nicer" when then making order monotone
 168                  * w.r.t. a simplicial function. */
 169                 void randomize_order() {
 170                     std::random_shuffle(order.begin(), order.end());
 171                     restore_revorder_from_order();
 172                 }
 173
 174                 /** Get the boundary matrix of the complex according to this
 175                  * order. Requires that is_permutation() gives true. */
 176                 boundary_matrix get_boundary_matrix() const {
 177                     assert(is_permutation());
 178                     boundary_matrix mat(size());
 179
 180                     for (unsigned c=0; c < size(); ++c)
 181                         for(unsigned r=0; r < get_simplex(c).face_count(); ++r)
 182                             mat.set(complex_to_order(get_simplex(c).faces[r]), c, true);
 183
 184                     return mat;
 185                 }
 186
 187                 /** Get the order of the simplices. */
 188                 const order_type& get_order() const {
 189                     return order;
 190                 }
 191
 192                 /** Set the order of the simplices. Requires that order has the
 193                  * same size as the underlying complex. Requires that the
 194                  * is_permutation() gives true after setting the order. */
 195                 void set_order(const order_type& o) {
 196                     assert(o.size() == c.size());
 197
 198                     order = o;
 199                     restore_revorder_from_order();
 200                 }
 201
 202             private:
 203                 /** Reconstruct 'revorder' by inverting the permutation given by 'order'. */
 204                 void restore_revorder_from_order() {
 205                     assert(is_permutation());
 206
 207                     // Make revorder * order the identity permutation
 208                     for (unsigned i=0; i < size(); ++i)
 209                         revorder[order[i]] = i;
 210                 }
 211
 212                 /** The complex of which we consider a simplex order. */
 213                 const scomplex &c;
 214
 215                 /** The i-th simplex in order is the order[i]-th simplex of the
 216                  * complex. 'order' can be seen as a permutation of the
 217                  * simplices saved in 'c'. */
 218                 order_type order;
 219
 220                 /** The i-th simplex in the complex is the revorder[i]-th
 221                  * simplex in order. 'revorder' can be seen as the inverse
 222                  * permutation saved in 'order'. */
 223                 order_type revorder;
 224         };
 225
 226     public:
 227         simplicial_complex() {
 228             // Add the one and only minus-one dimensional simplex
 229             simplices.push_back(simplex::create_minusonedim_simplex());
 230         }
 231
 232         /** Remove all simplices except the dummy simplex */
 233         void clear() {
 234             simplices.resize(1);
 235         }
 236
 237         /** Return number of simplices. */
 238         size_t size() const {
 239             return simplices.size();
 240         }
 241
 242         /** Add a simplex to the complex. The dimension of the faces must be
 243          * dim-1, and they must already be part of the complex. Returns the
 244          * index of the added simplex. */
 245         index_type add_simplex(int dim, index_type* faces) {
 246             return add_simplex(simplex(dim, faces));
 247         }
 248
 249         /** Add a simplex to the complex. The dimension of the faces must be
 250          * dim-1, and they must already be part of the complex. Returns the
 251          * index of the added simplex. The dimension must be at least zero. */
 252         index_type add_simplex(simplex s) {
 253             assert(s.dim >= 0);
 254
 255             // Check requirements for faces
 256             for (unsigned i=0; i < s.face_count(); ++i) {
 257                 // Faces are already in complex.
 258                 assert(s.faces[i] < size());
 259                 // Faces have dimension dim-1
 260                 assert(simplices[s.faces[i]].dim == s.dim-1);
 261             }
 262
 263             // index_type must be large enough
 264             assert(simplices.size() < std::numeric_limits<IT>::max());
 265
 266             index_type idx = simplices.size();
 267             simplices.push_back(s);
 268             return idx;
 269         }
 270
 271         /** Add an array of simplices */
 272         void add_simplices(simplex* sarray, size_t count) {
 273             for (unsigned i=0; i < count; ++i)
 274                 add_simplex(sarray[i]);
 275         }
 276
 277         /** Get idx-th simplex */
 278         const simplex& get_simplex(index_type idx) const {
 279             assert(0 <= idx && idx < size());
 280             return simplices.at(idx);
 281         }
 282
 283         /** Return true iff for each simplex i with dimension dim it holds that
 284          * the faces of i are contained in the complex and have dimension dim-1. */
 285         bool is_complex() const {
 286             for (unsigned i=0; i < size(); ++i) {
 287
 288                 const simplex &s = simplices[i];
 289                 for (unsigned f=0; f < s.face_count(); ++f) {
 290
 291                     if (s.faces[f] >= size())
 292                         return false;
 293
 294                     const simplex &face = simplices[s.faces[f]];
 295                     if (face.dim != s.dim-1)
 296                         return false;
 297                 }
 298             }
 299             return true;
 300         }
 301
 302     private:
 303         /** A list of simplices */
 304         std::vector<simplex> simplices;
 305 };
 306
 307 }
 308
 309
 310 #endif