Sanitize index check
[libstick.git] / include / libstick-0.1 / simplicialcomplex.h
index 707fd4f3240ae7fe71c4c90e7a4d1ab25926689e..1695e25abc6fbe611bb7fa9bfc27e12a618c5f03 100644 (file)
@@ -10,7 +10,7 @@
 
 #include <iostream>
 
-#include <libstick-0.1/booleanmatrix.h>
+#include "booleanmatrix.h"
 
 
 namespace libstick {
@@ -23,7 +23,7 @@ namespace libstick {
  * Each 0-dimensional simplex automatically has this simplex as its face.
  * Consequently, the innner class simplex_order gives the extended boundary
  * matrix. */
-template<int MAXDIM, class IT=uint32_t, class VT=double>
+template<int MAXDIM, class IT, class VT>
 class simplicial_complex {
 
     public:
@@ -36,7 +36,7 @@ class simplicial_complex {
         typedef VT value_type;
 
         /** A simplex of the complex. */
-        struct Simplex {
+        struct simplex {
             /** Dimension of the simplex. */
             int dim;
             /** The indices of the faces of the simplex. */
@@ -47,10 +47,10 @@ class simplicial_complex {
             /** Create a new simplex with dimension 'dim', (dim+1)-faces and
              * its value. If simpley is 0-dimensional, its face is
              * automatically set to one (-1)-dimensional simplex. */
-            static Simplex create(int dim, index_type* faces, value_type value) {
+            static simplex create(int dim, index_type* faces, value_type value) {
                 assert(0 <= dim && dim <= MAXDIM);
 
-                Simplex s;
+                simplex s;
                 s.dim = dim;
                 s.value = value;
 
@@ -63,8 +63,8 @@ class simplicial_complex {
             }
 
             /** Create a (-1)-dimensional simplex. It has the lowest possible value. */
-            static Simplex create_minusonedim_simplex() {
-                Simplex s;
+            static simplex create_minusonedim_simplex() {
+                simplex s;
 
                 s.dim = -1;
                 s.faces[0] = 0;
@@ -92,7 +92,7 @@ class simplicial_complex {
         class simplex_order {
 
             public:
-                typedef boolean_colrowmatrix<IT> boundary_matrix;
+                typedef boolean_colmatrix<IT> boundary_matrix;
 
                 /** Create a standard order of the complex c, i.e., the identity permutation. */
                 simplex_order(const simplcompltype &c) :
@@ -116,11 +116,15 @@ class simplicial_complex {
                 }
 
                 /** Get i-th simplex in the simplex order. */
-                const Simplex& get_simplex(size_t i) const {
-                    assert(i < size());
+                const simplex& get_simplex(index_type i) const {
+                    assert(0 <= i && i < size());
                     return c.simplices[order.at(i)];
                 }
 
+                const simplcompltype& get_complex() const {
+                    return c;
+                }
+
                 /** Returns true iff the faces of simplex i are before i in this order. */
                 bool is_filtration() const {
                     assert(size() == c.size());
@@ -145,13 +149,22 @@ class simplicial_complex {
                     return is_filtration();
                 }
 
+                /** Randomize order. It has hardly any impact on runtime, but
+                 * it makes cycles "nicer" when the simplice's function values
+                 * are constant.
+                 * */
+                void randomize_order() {
+                    std::random_shuffle(order.begin(), order.end());
+                    restore_revorder_from_order();
+                }
+
                 /** Sort simplices such that is_monotone() gives true. This
                  * requires that the complex's is_monotone() gave true
                  * beforehand.*/
                 void make_monotone_filtration() {
                     assert(c.is_monotone());
 
-                    sort(order.begin(), order.end(), cmp_monotone_filtration(c));
+                    std::sort(order.begin(), order.end(), cmp_monotone_filtration(c));
                     restore_revorder_from_order();
 
                     assert(c.is_monotone());
@@ -195,7 +208,12 @@ class simplicial_complex {
     public:
         simplicial_complex() {
             // Add the one minus-one dimensional simplex
-            add_simplex(Simplex::create_minusonedim_simplex());
+            add_simplex(simplex::create_minusonedim_simplex());
+        }
+
+        /** Remove all simplices except the dummy simplex */
+        void clear() {
+            simplices.resize(1);
         }
 
         /** Return number of simplices. */
@@ -204,14 +222,31 @@ class simplicial_complex {
         }
 
         /** Add a simplex to the complex. The dimension of the faces must be
-         * dim-1, and they must already be part of the complex. */
-        void add_simplex(int dim, index_type* faces, value_type value) {
-            add_simplex(Simplex::create(dim, faces, value));
+         * dim-1, and they must already be part of the complex. Returns the
+         * index of the added simplex. */
+        index_type add_simplex(int dim, index_type* faces, value_type value) {
+            return add_simplex(simplex::create(dim, faces, value));
+        }
+
+        /** Add a simplex to the complex of at least dimension 1. The dimension
+         * of the faces must be dim-1, and they must already be part of the
+         * complex. The value of the simplex is set to the maximum value of its
+         * faces. Returns the index of the added simplex. */
+        index_type add_simplex(int dim, index_type* faces) {
+            assert(dim >= 1);
+
+            // Get max value of its faces
+            VT value = simplices[faces[0]].value;
+            for (size_t i=0; i < simplex::face_count_bydim(dim); ++i)
+                value = std::max(value, simplices[faces[i]].value);
+
+            return add_simplex(dim, faces, value);
         }
 
         /** Add a simplex to the complex. The dimension of the faces must be
-         * dim-1, and they must already be part of the complex. */
-        void add_simplex(Simplex s) {
+         * dim-1, and they must already be part of the complex. Returns the
+         * index of the added simplex. */
+        index_type add_simplex(simplex s) {
             // Check requirements for faces
             for (unsigned i=0; i < s.face_count(); ++i) {
                 // Faces are already in complex.
@@ -220,11 +255,16 @@ class simplicial_complex {
                 assert(simplices[s.faces[i]].dim == s.dim-1);
             }
 
+            // index_type must be large enough
+            assert(simplices.size() < std::numeric_limits<IT>::max());
+
+            index_type idx = simplices.size();
             simplices.push_back(s);
+            return idx;
         }
 
         /** Add an array of simplices */
-        void add_simplices(Simplex* sarray, size_t count) {
+        void add_simplices(simplex* sarray, size_t count) {
             for (unsigned i=0; i < count; ++i)
                 add_simplex(sarray[i]);
         }
@@ -234,13 +274,13 @@ class simplicial_complex {
         bool is_complex() const {
             for (unsigned i=0; i < size(); ++i) {
 
-                const Simplex &s = simplices[i];
+                const simplex &s = simplices[i];
                 for (unsigned f=0; f < s.face_count(); ++f) {
 
                     if (s.faces[f] >= size())
                         return false;
 
-                    const Simplex &face = simplices[s.faces[f]];
+                    const simplex &face = simplices[s.faces[f]];
                     if (face.dim != s.dim-1)
                         return false;
                 }
@@ -254,7 +294,7 @@ class simplicial_complex {
         bool is_monotone() const {
             assert(is_complex());
 
-            typename std::vector<Simplex>::const_iterator it = ++simplices.begin();
+            typename std::vector<simplex>::const_iterator it = ++simplices.begin();
             for (; it != simplices.end(); ++it)
                 for (unsigned f=0; f < it->face_count(); ++f)
                     if (simplices[it->faces[f]].value > it->value)
@@ -275,8 +315,8 @@ class simplicial_complex {
                 }
 
             bool operator()(index_type i, index_type j) {
-                const Simplex& si = c.simplices[i];
-                const Simplex& sj = c.simplices[j];
+                const simplex& si = c.simplices[i];
+                const simplex& sj = c.simplices[j];
 
                 if (si.value < sj.value)
                     return true;
@@ -289,7 +329,7 @@ class simplicial_complex {
 
     public:
         /** A list of simplices */
-        std::vector<Simplex> simplices;
+        std::vector<simplex> simplices;
 };
 
 }