Add matrix reduction code

[libstick.git] / include / libstick-0.1 / simplicialcomplex.h

diff --git a/include/libstick-0.1/simplicialcomplex.h b/include/libstick-0.1/simplicialcomplex.h
index d829721c02da3515a354ec89f034933ecb90e46a..07080ed76a23617e2c87b99ca7aa3dc998640b2f 100644 (file)
--- a/include/libstick-0.1/simplicialcomplex.h
+++ b/include/libstick-0.1/simplicialcomplex.h
@@ -6,6 +6,7 @@
 #include <cstring>
 #include <algorithm>
 #include <vector>
+#include <limits>
 
 #include <iostream>
 
@@ -17,8 +18,12 @@ namespace libstick {
 /** A simplicial complex is a std::vector of simplices such that each face is
  * also part of the complex. Every simplex has dimension at most MAXDIM. The
  * indices of simplices resp. their faces are of type IT. To each simplex a
- * value is assigend, which is of type VT. */
-template<size_t MAXDIM, class IT=uint32_t, class VT=double>
+ * value is assigend, which is of type VT. When a SimplicialComplex is
+ * instantiated, a single (-1) dimensional simplex is automatically created.
+ * Each 0-dimensional simplex automatically has this simplex as its face.
+ * Consequently, the innner class SimplexOrder gives the extended boundary
+ * matrix. */
+template<int MAXDIM, class IT=uint32_t, class VT=double>
 class SimplicialComplex {
 
     public:
@@ -33,34 +38,51 @@ class SimplicialComplex {
         /** A simplex of the complex. */
         struct Simplex {
             /** Dimension of the simplex. */
-            size_t dim;
+            int dim;
             /** The indices of the faces of the simplex. */
             index_type faces[MAXDIM+1];
             /** The value of the simplex. */
             value_type value;
 
             /** Create a new simplex with dimension 'dim', (dim+1)-faces and
-             * its value. */
-            static Simplex create(size_t dim, index_type* faces, value_type value) {
-                assert(dim <= MAXDIM);
+             * 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) {
+                assert(0 <= dim && dim <= MAXDIM);
 
                 Simplex s;
                 s.dim = dim;
                 s.value = value;
-                memcpy(s.faces, faces, face_count_bydim(dim)*sizeof(index_type));
+
+                if (dim > 0)
+                    memcpy(s.faces, faces, faceCountByDim(dim)*sizeof(index_type));
+                else
+                    s.faces[0] = 0;
+
+                return s;
+            }
+
+            /** Create a (-1)-dimensional simplex. It has the lowest possible value. */
+            static Simplex create_minusonedim_simplex() {
+                Simplex s;
+
+                s.dim = -1;
+                s.faces[0] = 0;
+                s.value = std::numeric_limits<VT>::has_infinity
+                            ? -std::numeric_limits<VT>::infinity()
+                            : std::numeric_limits<VT>::min();
 
                 return s;
             }
 
             /** Get number of faces. */
-            size_t face_count() const {
-                return face_count_bydim(dim);
+            size_t faceCount() const {
+                return faceCountByDim(dim);
             }
 
             /** Get number of faces of a dim-dimensional simplex. */
-            static size_t face_count_bydim(size_t dim) {
-                if (dim == 0)
-                    return 0;
+            static size_t faceCountByDim(int dim) {
+                assert(-1 <= dim && dim <= MAXDIM);
                 return dim + 1;
             }
         };
@@ -104,7 +126,7 @@ class SimplicialComplex {
                     assert(size() == c.size());
 
                     for (unsigned i=0; i < size(); ++i)
-                        for (unsigned f=0; f < getSimplex(i).face_count(); ++f)
+                        for (unsigned f=0; f < getSimplex(i).faceCount(); ++f)
                             if (revorder[getSimplex(i).faces[f]] >= i)
                                 return false;
 
@@ -142,7 +164,7 @@ class SimplicialComplex {
                     BoundaryMatrix mat(size());
 
                     for (unsigned c=0; c < size(); ++c)
-                        for(unsigned r=0; r < getSimplex(c).face_count(); ++r)
+                        for(unsigned r=0; r < getSimplex(c).faceCount(); ++r)
                             mat.set(revorder[getSimplex(c).faces[r]], c, true);
 
                     return mat;
@@ -171,6 +193,11 @@ class SimplicialComplex {
         };
 
     public:
+        SimplicialComplex() {
+            // Add the one minus-one dimensional simplex
+            addSimplex(Simplex::create_minusonedim_simplex());
+        }
+
         /** Return number of simplices. */
         size_t size() const {
             return simplices.size();
@@ -178,7 +205,7 @@ class SimplicialComplex {
 
         /** 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 addSimplex(size_t dim, index_type* faces, value_type value) {
+        void addSimplex(int dim, index_type* faces, value_type value) {
             addSimplex(Simplex::create(dim, faces, value));
         }
 
@@ -186,7 +213,7 @@ class SimplicialComplex {
          * dim-1, and they must already be part of the complex. */
         void addSimplex(Simplex s) {
             // Check requirements for faces
-            for (unsigned i=0; i < s.face_count(); ++i) {
+            for (unsigned i=0; i < s.faceCount(); ++i) {
                 // Faces are already in complex.
                 assert(s.faces[i] < size());
                 // Faces have dimension dim-1
@@ -199,11 +226,10 @@ class SimplicialComplex {
         /** Return true iff for each simplex i with dimension dim it holds that
          * the faces of i are contained in the complex and have dimension dim-1. */
         bool isComplex() const {
-            typename std::vector<Simplex>::const_iterator it = ++simplices.begin();
             for (unsigned i=0; i < size(); ++i) {
 
                 const Simplex &s = simplices[i];
-                for (unsigned f=0; f < s.face_count(); ++f) {
+                for (unsigned f=0; f < s.faceCount(); ++f) {
 
                     if (s.faces[f] >= size())
                         return false;
@@ -224,7 +250,7 @@ class SimplicialComplex {
 
             typename std::vector<Simplex>::const_iterator it = ++simplices.begin();
             for (; it != simplices.end(); ++it)
-                for (unsigned f=0; f < it->face_count(); ++f)
+                for (unsigned f=0; f < it->faceCount(); ++f)
                     if (simplices[it->faces[f]].value > it->value)
                         return false;