include/libstick-0.1/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 std::vector 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 resp. their faces are of type IT. To each simplex a
  21  * value is assigend, which is of type VT. When a simplicial_complex is
  22  * instantiated, a single (-1) dimensional simplex is automatically created.
  23  * Each 0-dimensional simplex automatically has this simplex as its face.
  24  * Consequently, the innner class simplex_order gives the extended boundary
  25  * matrix. */
  26 template<int MAXDIM, class IT, class VT>
  27 class simplicial_complex {
  28
  29     public:
  30         /** The type of this class. */
  31         typedef simplicial_complex<MAXDIM, IT, VT> simplcompltype;
  32         /** Type of indices of simplices. */
  33         typedef IT index_type;
  34         /** To every simplex a function value is assigned according to which a
  35          * filtration is considered. This is the value type of the function. */
  36         typedef VT value_type;
  37
  38         /** A simplex of the complex. */
  39         struct simplex {
  40             /** Dimension of the simplex. */
  41             int dim;
  42             /** The indices of the faces of the simplex. */
  43             index_type faces[MAXDIM+1];
  44             /** The value of the simplex. */
  45             value_type value;
  46
  47             /** Create a new simplex with dimension 'dim', (dim+1)-faces and
  48              * its value. If simpley is 0-dimensional, its face is
  49              * automatically set to one (-1)-dimensional simplex. */
  50             static simplex create(int dim, index_type* faces, value_type value) {
  51                 assert(0 <= dim && dim <= MAXDIM);
  52
  53                 simplex s;
  54                 s.dim = dim;
  55                 s.value = value;
  56
  57                 if (dim > 0)
  58                     memcpy(s.faces, faces, face_count_bydim(dim)*sizeof(index_type));
  59                 else
  60                     s.faces[0] = 0;
  61
  62                 return s;
  63             }
  64
  65             /** Create a (-1)-dimensional simplex. It has the lowest possible value. */
  66             static simplex create_minusonedim_simplex() {
  67                 simplex s;
  68
  69                 s.dim = -1;
  70                 s.faces[0] = 0;
  71                 s.value = std::numeric_limits<VT>::has_infinity
  72                             ? -std::numeric_limits<VT>::infinity()
  73                             : std::numeric_limits<VT>::min();
  74
  75                 return s;
  76             }
  77
  78             /** Get number of faces. */
  79             size_t face_count() const {
  80                 return face_count_bydim(dim);
  81             }
  82
  83             /** Get number of faces of a dim-dimensional simplex. */
  84             static size_t face_count_bydim(int dim) {
  85                 assert(-1 <= dim && dim <= MAXDIM);
  86                 return dim + 1;
  87             }
  88         };
  89
  90         /** An order of the simplices of complex c. An order can be interpreted
  91          * as a permuation of the complex's std::vector of simplices.  */
  92         class simplex_order {
  93
  94             public:
  95                 typedef boolean_colmatrix<IT> boundary_matrix;
  96
  97                 /** Create a standard order of the complex c, i.e., the identity permutation. */
  98                 simplex_order(const simplcompltype &c) :
  99                     c(c)
 100                     {
 101                     reset();
 102                 }
 103
 104                 /** Reset order to the identity permutation of the complex's simplices. */
 105                 void reset() {
 106                     order.clear();
 107                     for (unsigned i=0; i < c.size(); ++i)
 108                         order.push_back(i);
 109                     revorder = order;
 110                 }
 111
 112                 /** Return number of simplices. */
 113                 size_t size() const {
 114                     assert(order.size() == revorder.size());
 115                     return order.size();
 116                 }
 117
 118                 /** Get i-th simplex in the simplex order. */
 119                 const simplex& get_simplex(index_type i) const {
 120                     assert(0 <= i && i < size());
 121                     return c.simplices[order.at(i)];
 122                 }
 123
 124                 const simplcompltype& get_complex() const {
 125                     return c;
 126                 }
 127
 128                 /** Returns true iff the faces of simplex i are before i in this order. */
 129                 bool is_filtration() const {
 130                     assert(size() == c.size());
 131
 132                     for (unsigned i=0; i < size(); ++i)
 133                         for (unsigned f=0; f < get_simplex(i).face_count(); ++f)
 134                             if (revorder[get_simplex(i).faces[f]] >= i)
 135                                 return false;
 136
 137                     return true;
 138                 }
 139
 140                 /** Returns true iff is_filtration() gives true and values of simplices
 141                  * are monotone w.r.t. this order of simplices. */
 142                 bool is_monotone() const {
 143                     assert(size() == c.size());
 144
 145                     for (unsigned i=1; i < size(); ++i)
 146                         if (get_simplex(i-1).value > get_simplex(i).value)
 147                             return false;
 148
 149                     return is_filtration();
 150                 }
 151
 152                 /** Randomize order. It has hardly any impact on runtime, but
 153                  * it makes cycles "nicer" when the simplice's function values
 154                  * are constant.
 155                  * */
 156                 void randomize_order() {
 157                     std::random_shuffle(order.begin(), order.end());
 158                     restore_revorder_from_order();
 159                 }
 160
 161                 /** Sort simplices such that is_monotone() gives true. This
 162                  * requires that the complex's is_monotone() gave true
 163                  * beforehand.*/
 164                 void make_monotone_filtration() {
 165                     assert(c.is_monotone());
 166
 167                     std::sort(order.begin(), order.end(), cmp_monotone_filtration(c));
 168                     restore_revorder_from_order();
 169
 170                     assert(c.is_monotone());
 171                     assert(is_filtration());
 172                     assert(is_monotone());
 173                 }
 174
 175                 /** Get the boundary matrix of the complex according to this order. */
 176                 boundary_matrix get_boundary_matrix() const {
 177                     boundary_matrix mat(size());
 178
 179                     for (unsigned c=0; c < size(); ++c)
 180                         for(unsigned r=0; r < get_simplex(c).face_count(); ++r)
 181                             mat.set(revorder[get_simplex(c).faces[r]], c, true);
 182
 183                     return mat;
 184                 }
 185
 186             private:
 187                 /** Reconstruct 'revorder' by inverting the permutation given by 'order'. */
 188                 void restore_revorder_from_order() {
 189                     // Make revorder * order the identity permutation
 190                     for (unsigned i=0; i < size(); ++i)
 191                         revorder[order[i]] = i;
 192                 }
 193
 194                 /** The complex of which we consider a simplex order. */
 195                 const simplcompltype &c;
 196
 197                 /** The i-th simplex in order is the order[i]-th simplex of the
 198                  * complex. 'order' can be seen as a permutation of the
 199                  * simplices saved in 'c'. */
 200                 std::vector<index_type> order;
 201
 202                 /** The i-th simplex in the complex is the revorder[i]-th
 203                  * simplex in order. 'revorder' can be seen as the inverse
 204                  * permutation saved in 'order'. */
 205                 std::vector<index_type> revorder;
 206         };
 207
 208     public:
 209         simplicial_complex() {
 210             // Add the one minus-one dimensional simplex
 211             add_simplex(simplex::create_minusonedim_simplex());
 212         }
 213
 214         /** Remove all simplices except the dummy simplex */
 215         void clear() {
 216             simplices.resize(1);
 217         }
 218
 219         /** Return number of simplices. */
 220         size_t size() const {
 221             return simplices.size();
 222         }
 223
 224         /** Add a simplex to the complex. The dimension of the faces must be
 225          * dim-1, and they must already be part of the complex. Returns the
 226          * index of the added simplex. */
 227         index_type add_simplex(int dim, index_type* faces, value_type value) {
 228             return add_simplex(simplex::create(dim, faces, value));
 229         }
 230
 231         /** Add a simplex to the complex of at least dimension 1. The dimension
 232          * of the faces must be dim-1, and they must already be part of the
 233          * complex. The value of the simplex is set to the maximum value of its
 234          * faces. Returns the index of the added simplex. */
 235         index_type add_simplex(int dim, index_type* faces) {
 236             assert(dim >= 1);
 237
 238             // Get max value of its faces
 239             VT value = simplices[faces[0]].value;
 240             for (size_t i=0; i < simplex::face_count_bydim(dim); ++i)
 241                 value = std::max(value, simplices[faces[i]].value);
 242
 243             return add_simplex(dim, faces, value);
 244         }
 245
 246         /** Add a simplex to the complex. The dimension of the faces must be
 247          * dim-1, and they must already be part of the complex. Returns the
 248          * index of the added simplex. */
 249         index_type add_simplex(simplex s) {
 250             // Check requirements for faces
 251             for (unsigned i=0; i < s.face_count(); ++i) {
 252                 // Faces are already in complex.
 253                 assert(s.faces[i] < size());
 254                 // Faces have dimension dim-1
 255                 assert(simplices[s.faces[i]].dim == s.dim-1);
 256             }
 257
 258             // index_type must be large enough
 259             assert(simplices.size() < std::numeric_limits<IT>::max());
 260
 261             index_type idx = simplices.size();
 262             simplices.push_back(s);
 263             return idx;
 264         }
 265
 266         /** Add an array of simplices */
 267         void add_simplices(simplex* sarray, size_t count) {
 268             for (unsigned i=0; i < count; ++i)
 269                 add_simplex(sarray[i]);
 270         }
 271
 272         /** Return true iff for each simplex i with dimension dim it holds that
 273          * the faces of i are contained in the complex and have dimension dim-1. */
 274         bool is_complex() const {
 275             for (unsigned i=0; i < size(); ++i) {
 276
 277                 const simplex &s = simplices[i];
 278                 for (unsigned f=0; f < s.face_count(); ++f) {
 279
 280                     if (s.faces[f] >= size())
 281                         return false;
 282
 283                     const simplex &face = simplices[s.faces[f]];
 284                     if (face.dim != s.dim-1)
 285                         return false;
 286                 }
 287             }
 288             return true;
 289         }
 290
 291         /** Returns true iff simplex's values are monotone w.r.t.
 292          * face-inclusion, i.e., for each simplex its value is not smaller than
 293          * the values of its faces. Requires that is_complex() gives true. */
 294         bool is_monotone() const {
 295             assert(is_complex());
 296
 297             typename std::vector<simplex>::const_iterator it = ++simplices.begin();
 298             for (; it != simplices.end(); ++it)
 299                 for (unsigned f=0; f < it->face_count(); ++f)
 300                     if (simplices[it->faces[f]].value > it->value)
 301                         return false;
 302
 303             return true;
 304         }
 305
 306     private:
 307         /** Compares (operator<) two simplices (i.e. indices) in a
 308          * simplex_order w.r.t. lexicographical order on (value,
 309          * dimension)-tuples. */
 310         struct cmp_monotone_filtration {
 311             const simplicial_complex &c;
 312
 313             cmp_monotone_filtration(const simplicial_complex &c) :
 314                 c(c){
 315                 }
 316
 317             bool operator()(index_type i, index_type j) {
 318                 const simplex& si = c.simplices[i];
 319                 const simplex& sj = c.simplices[j];
 320
 321                 if (si.value < sj.value)
 322                     return true;
 323                 else if (si.value == sj.value)
 324                     return si.dim < sj.dim;
 325                 else
 326                     return false;
 327             }
 328         };
 329
 330     public:
 331         /** A list of simplices */
 332         std::vector<simplex> simplices;
 333 };
 334
 335 }
 336
 337
 338 #endif