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
  10 #include <iostream>
  11
  12 #include <libstick-0.1/booleanmatrix.h>
  13
  14
  15 namespace libstick {
  16
  17 /** A simplicial complex is a std::vector of simplices such that each face is
  18  * also part of the complex. Every simplex has dimension at most MAXDIM. The
  19  * indices of simplices resp. their faces are of type IT. To each simplex a
  20  * value is assigend, which is of type VT. */
  21 template<size_t MAXDIM, class IT=uint32_t, class VT=double>
  22 class SimplicialComplex {
  23
  24     public:
  25         /** The type of this class. */
  26         typedef SimplicialComplex<MAXDIM, IT, VT> simplcompltype;
  27         /** Type of indices of simplices. */
  28         typedef IT index_type;
  29         /** To every simplex a function value is assigned according to which a
  30          * filtration is considered. This is the value type of the function. */
  31         typedef VT value_type;
  32
  33         /** A simplex of the complex. */
  34         struct Simplex {
  35             /** Dimension of the simplex. */
  36             size_t dim;
  37             /** The indices of the faces of the simplex. */
  38             index_type faces[MAXDIM+1];
  39             /** The value of the simplex. */
  40             value_type value;
  41
  42             /** Create a new simplex with dimension 'dim', (dim+1)-faces and
  43              * its value. */
  44             static Simplex create(size_t dim, index_type* faces, value_type value) {
  45                 assert(dim <= MAXDIM);
  46
  47                 Simplex s;
  48                 s.dim = dim;
  49                 s.value = value;
  50                 memcpy(s.faces, faces, face_count_bydim(dim)*sizeof(index_type));
  51
  52                 return s;
  53             }
  54
  55             /** Get number of faces. */
  56             size_t face_count() const {
  57                 return face_count_bydim(dim);
  58             }
  59
  60             /** Get number of faces of a dim-dimensional simplex. */
  61             static size_t face_count_bydim(size_t dim) {
  62                 if (dim == 0)
  63                     return 0;
  64                 return dim + 1;
  65             }
  66         };
  67
  68         /** An order of the simplices of complex c. An order can be interpreted
  69          * as a permuation of the complex's std::vector of simplices.  */
  70         class SimplexOrder {
  71
  72             public:
  73                 typedef BooleanColRowMatrix<IT> BoundaryMatrix;
  74
  75                 /** Create a standard order of the complex c, i.e., the identity permutation. */
  76                 SimplexOrder(const simplcompltype &c) :
  77                     c(c)
  78                     {
  79                     reset();
  80                 }
  81
  82                 /** Reset order to the identity permutation of the complex's simplices. */
  83                 void reset() {
  84                     order.clear();
  85                     for (unsigned i=0; i < c.size(); ++i)
  86                         order.push_back(i);
  87                     revorder = order;
  88                 }
  89
  90                 /** Return number of simplices. */
  91                 size_t size() const {
  92                     assert(order.size() == revorder.size());
  93                     return order.size();
  94                 }
  95
  96                 /** Get i-th simplex in the simplex order. */
  97                 const Simplex& getSimplex(size_t i) const {
  98                     assert(i < size());
  99                     return c.simplices[order.at(i)];
 100                 }
 101
 102                 /** Returns true iff the faces of simplex i are before i in this order. */
 103                 bool isFiltration() const {
 104                     assert(size() == c.size());
 105
 106                     for (unsigned i=0; i < size(); ++i)
 107                         for (unsigned f=0; f < getSimplex(i).face_count(); ++f)
 108                             if (revorder[getSimplex(i).faces[f]] >= i)
 109                                 return false;
 110
 111                     return true;
 112                 }
 113
 114                 /** Returns true iff isFiltration() gives true and values of simplices
 115                  * are monotone w.r.t. this order of simplices. */
 116                 bool isMonotone() const {
 117                     assert(size() == c.size());
 118
 119                     for (unsigned i=1; i < size(); ++i)
 120                         if (getSimplex(i-1).value > getSimplex(i).value)
 121                             return false;
 122
 123                     return isFiltration();
 124                 }
 125
 126                 /** Sort simplices such that isMonotone() gives true. This
 127                  * requires that the complex's isMonotone() gave true
 128                  * beforehand.*/
 129                 void makeMonotoneFiltration() {
 130                     assert(c.isMonotone());
 131
 132                     sort(order.begin(), order.end(), cmpMonotoneFiltration(c));
 133                     restoreRevorderFromOrder();
 134
 135                     assert(c.isMonotone());
 136                     assert(isFiltration());
 137                     assert(isMonotone());
 138                 }
 139
 140                 /** Get the boundary matrix of the complex according to this order. */
 141                 BoundaryMatrix getBoundaryMatrix() const {
 142                     BoundaryMatrix mat(size());
 143
 144                     for (unsigned c=0; c < size(); ++c)
 145                         for(unsigned r=0; r < getSimplex(c).face_count(); ++r)
 146                             mat.set(revorder[getSimplex(c).faces[r]], c, true);
 147
 148                     return mat;
 149                 }
 150
 151             private:
 152                 /** Reconstruct 'revorder' by inverting the permutation given by 'order'. */
 153                 void restoreRevorderFromOrder() {
 154                     // Make revorder * order the identity permutation
 155                     for (unsigned i=0; i < size(); ++i)
 156                         revorder[order[i]] = i;
 157                 }
 158
 159                 /** The complex of which we consider a simplex order. */
 160                 const simplcompltype &c;
 161
 162                 /** The i-th simplex in order is the order[i]-th simplex of the
 163                  * complex. 'order' can be seen as a permutation of the
 164                  * simplices saved in 'c'. */
 165                 std::vector<index_type> order;
 166
 167                 /** The i-th simplex in the complex is the revorder[i]-th
 168                  * simplex in order. 'revorder' can be seen as the inverse
 169                  * permutation saved in 'order'. */
 170                 std::vector<index_type> revorder;
 171         };
 172
 173     public:
 174         /** Return number of simplices. */
 175         size_t size() const {
 176             return simplices.size();
 177         }
 178
 179         /** Add a simplex to the complex. The dimension of the faces must be
 180          * dim-1, and they must already be part of the complex. */
 181         void addSimplex(size_t dim, index_type* faces, value_type value) {
 182             addSimplex(Simplex::create(dim, faces, value));
 183         }
 184
 185         /** Add a simplex to the complex. The dimension of the faces must be
 186          * dim-1, and they must already be part of the complex. */
 187         void addSimplex(Simplex s) {
 188             // Check requirements for faces
 189             for (unsigned i=0; i < s.face_count(); ++i) {
 190                 // Faces are already in complex.
 191                 assert(s.faces[i] < size());
 192                 // Faces have dimension dim-1
 193                 assert(simplices[s.faces[i]].dim == s.dim-1);
 194             }
 195
 196             simplices.push_back(s);
 197         }
 198
 199         /** Return true iff for each simplex i with dimension dim it holds that
 200          * the faces of i are contained in the complex and have dimension dim-1. */
 201         bool isComplex() const {
 202             typename std::vector<Simplex>::const_iterator it = ++simplices.begin();
 203             for (unsigned i=0; i < size(); ++i) {
 204
 205                 const Simplex &s = simplices[i];
 206                 for (unsigned f=0; f < s.face_count(); ++f) {
 207
 208                     if (s.faces[f] >= size())
 209                         return false;
 210
 211                     const Simplex &face = simplices[s.faces[f]];
 212                     if (face.dim != s.dim-1)
 213                         return false;
 214                 }
 215             }
 216             return true;
 217         }
 218
 219         /** Returns true iff simplex's values are monotone w.r.t.
 220          * face-inclusion, i.e., for each simplex its value is not smaller than
 221          * the values of its faces. Requires that isComplex() gives true. */
 222         bool isMonotone() const {
 223             assert(isComplex());
 224
 225             typename std::vector<Simplex>::const_iterator it = ++simplices.begin();
 226             for (; it != simplices.end(); ++it)
 227                 for (unsigned f=0; f < it->face_count(); ++f)
 228                     if (simplices[it->faces[f]].value > it->value)
 229                         return false;
 230
 231             return true;
 232         }
 233
 234     private:
 235         /** Compares (operator<) two simplices (i.e. indices) in a
 236          * SimplexOrder w.r.t. lexicographical order on (value,
 237          * dimension)-tuples. */
 238         struct cmpMonotoneFiltration {
 239             const SimplicialComplex &c;
 240
 241             cmpMonotoneFiltration(const SimplicialComplex &c) :
 242                 c(c){
 243                 }
 244
 245             bool operator()(index_type i, index_type j) {
 246                 const Simplex& si = c.simplices[i];
 247                 const Simplex& sj = c.simplices[j];
 248
 249                 if (si.value < sj.value)
 250                     return true;
 251                 else if (si.value == sj.value)
 252                     return si.dim < sj.dim;
 253                 else
 254                     return false;
 255             }
 256         };
 257
 258     public:
 259         /** A list of simplices */
 260         std::vector<Simplex> simplices;
 261 };
 262
 263 }
 264
 265
 266 #endif