X-Git-Url: https://git.sthu.org/?p=libstick.git;a=blobdiff_plain;f=include%2Flibstick%2Fsimplicialcomplex.h;h=6dd99a376f50b9392306941ed67524e7ad60a486;hp=1695e25abc6fbe611bb7fa9bfc27e12a618c5f03;hb=7186ae9cc97393496d7f0e0cb281524276f1e47c;hpb=d6673c66ae1cb2b421b6e0c4b6a2f6ea6f0fad98 diff --git a/include/libstick/simplicialcomplex.h b/include/libstick/simplicialcomplex.h index 1695e25..6dd99a3 100644 --- a/include/libstick/simplicialcomplex.h +++ b/include/libstick/simplicialcomplex.h @@ -15,63 +15,60 @@ namespace libstick { -/** A simplicial complex is a std::vector of simplices such that each face is +/** A simplicial complex is a collection 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. When a simplicial_complex 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 simplex_order gives the extended boundary - * matrix. */ -template + * indices of simplices and their faces are of type IT. When a + * simplicial_complex is instantiated, a single (-1) dimensional simplex is + * automatically created. (When we compute persistent homology, we actually + * compute reduced homology.) Each 0-dimensional simplex automatically has this + * simplex as its face. + * + * Based on a simplicial complex, we define a simplicial order. An order + * can be seen as a permutation of the simplices of the complex. If the order + * has the property that every prefix of the permutation is again a complex, + * i.e., all faces of all simplices of each prefix are contained in the prefix, + * then we call it a filtration. The innner class simplex_order gives the + * extended boundary matrix, which can then be used in the persistence class to + * compute persistent homology. + * + * One way to obtain a filtration is to define a monotone simplicial function + * on the complex. That is, each simplex gets a value of type VT assigned. Then + * one can obtain a sub-level set filtration from the complex w.r.t. to this + * simplicial function. + */ +template class simplicial_complex { public: /** The type of this class. */ - typedef simplicial_complex simplcompltype; + typedef simplicial_complex simplcompltype; /** Type of indices of simplices. */ typedef IT index_type; - /** To every simplex a function value is assigned according to which a - * filtration is considered. This is the value type of the function. */ - typedef VT value_type; - /** A simplex of the complex. */ + /** A simplex of the complex. This struct is an aggregate. */ struct simplex { /** Dimension of the simplex. */ 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. 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) { + /** Create a new simplex with dimension 'dim' and (dim+1)-faces. If + * the simplex is 0-dimensional, its face is automatically set to + * one (-1)-dimensional simplex. */ + static simplex create_simplex(int dim, index_type* faces) { assert(0 <= dim && dim <= MAXDIM); - simplex s; - s.dim = dim; - s.value = value; + simplex s = {dim, {0}}; - if (dim > 0) + if (dim >= 1) memcpy(s.faces, faces, face_count_bydim(dim)*sizeof(index_type)); - else - s.faces[0] = 0; return s; } - /** Create a (-1)-dimensional simplex. It has the lowest possible value. */ + /** Create a (-1)-dimensional simplex. */ static simplex create_minusonedim_simplex() { - simplex s; - - s.dim = -1; - s.faces[0] = 0; - s.value = std::numeric_limits::has_infinity - ? -std::numeric_limits::infinity() - : std::numeric_limits::min(); - + simplex s = {-1, {0}}; return s; } @@ -93,11 +90,11 @@ class simplicial_complex { public: typedef boolean_colmatrix boundary_matrix; + typedef std::vector order_type; /** Create a standard order of the complex c, i.e., the identity permutation. */ simplex_order(const simplcompltype &c) : - c(c) - { + c(c) { reset(); } @@ -117,75 +114,96 @@ class simplicial_complex { /** Get i-th simplex in the simplex order. */ const simplex& get_simplex(index_type i) const { + return c.simplices.at(order_to_complex(i)); + } + + /** Translate an index in the order to a simplex in the complex. */ + index_type order_to_complex(index_type i) const { + assert(0 <= i && i < size()); + return order.at(i); + } + + /** Translate an index in the complex to a simplex in the order. */ + index_type complex_to_order(index_type i) const { assert(0 <= i && i < size()); - return c.simplices[order.at(i)]; + return revorder.at(i); } + /** Return the underlying complex. */ 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()); + /** Check if this order is a permutation of the underlying + * complex. */ + bool is_permutation() const { + if (size() != c.size()) + return false; + // Sort the index array and check if every index is present. + order_type o = order; + std::sort(o.begin(), o.end()); for (unsigned i=0; i < size(); ++i) - for (unsigned f=0; f < get_simplex(i).face_count(); ++f) - if (revorder[get_simplex(i).faces[f]] >= i) - return false; + if (o.at(i) != i) + return false; return true; } - /** Returns true iff is_filtration() gives true and values of simplices - * are monotone w.r.t. this order of simplices. */ - bool is_monotone() const { + /** Returns true iff the faces of simplex i are before i in + * this order. */ + bool is_filtration() const { assert(size() == c.size()); - for (unsigned i=1; i < size(); ++i) - if (get_simplex(i-1).value > get_simplex(i).value) - return false; + for (unsigned i=0; i < size(); ++i) + for (unsigned f=0; f < get_simplex(i).face_count(); ++f) + if (complex_to_order(get_simplex(i).faces[f]) >= i) + return false; - return is_filtration(); + return true; } /** Randomize order. It has hardly any impact on runtime, but - * it makes cycles "nicer" when the simplice's function values - * are constant. - * */ + * it makes cycles "nicer" when then making order monotone + * w.r.t. a simplicial function. */ 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()); - - std::sort(order.begin(), order.end(), cmp_monotone_filtration(c)); - restore_revorder_from_order(); - - assert(c.is_monotone()); - assert(is_filtration()); - assert(is_monotone()); - } - - /** Get the boundary matrix of the complex according to this order. */ + /** Get the boundary matrix of the complex according to this + * order. Requires that is_permutation() gives true. */ boundary_matrix get_boundary_matrix() const { + assert(is_permutation()); boundary_matrix mat(size()); for (unsigned c=0; c < size(); ++c) for(unsigned r=0; r < get_simplex(c).face_count(); ++r) - mat.set(revorder[get_simplex(c).faces[r]], c, true); + mat.set(complex_to_order(get_simplex(c).faces[r]), c, true); return mat; } + /** Get the order of the simplices. */ + const order_type& get_order() const { + return order; + } + + /** Set the order of the simplices. Requires that order has the + * same size as the underlying complex. Requires that the + * is_permutation() gives true after setting the order. */ + void set_order(const order_type& o) { + assert(o.size() == c.size()); + + order = o; + restore_revorder_from_order(); + } + private: /** Reconstruct 'revorder' by inverting the permutation given by 'order'. */ void restore_revorder_from_order() { + assert(is_permutation()); + // Make revorder * order the identity permutation for (unsigned i=0; i < size(); ++i) revorder[order[i]] = i; @@ -197,18 +215,18 @@ class simplicial_complex { /** The i-th simplex in order is the order[i]-th simplex of the * complex. 'order' can be seen as a permutation of the * simplices saved in 'c'. */ - std::vector order; + order_type order; /** The i-th simplex in the complex is the revorder[i]-th * simplex in order. 'revorder' can be seen as the inverse * permutation saved in 'order'. */ - std::vector revorder; + order_type revorder; }; public: simplicial_complex() { - // Add the one minus-one dimensional simplex - add_simplex(simplex::create_minusonedim_simplex()); + // Add the one and only minus-one dimensional simplex + simplices.push_back(simplex::create_minusonedim_simplex()); } /** Remove all simplices except the dummy simplex */ @@ -224,29 +242,16 @@ 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. 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); + return add_simplex(simplex(dim, faces)); } /** Add a simplex to the complex. The dimension of the faces must be * dim-1, and they must already be part of the complex. Returns the - * index of the added simplex. */ + * index of the added simplex. The dimension must be at least zero. */ index_type add_simplex(simplex s) { + assert(s.dim >= 0); + // Check requirements for faces for (unsigned i=0; i < s.face_count(); ++i) { // Faces are already in complex. @@ -269,6 +274,12 @@ class simplicial_complex { add_simplex(sarray[i]); } + /** Get idx-th simplex */ + const simplex& get_simplex(index_type idx) const { + assert(0 <= idx && idx < size()); + return simplices.at(idx); + } + /** 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 is_complex() const { @@ -288,46 +299,7 @@ class simplicial_complex { return true; } - /** Returns true iff simplex's values are monotone w.r.t. - * face-inclusion, i.e., for each simplex its value is not smaller than - * the values of its faces. Requires that is_complex() gives true. */ - bool is_monotone() const { - assert(is_complex()); - - typename std::vector::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) - return false; - - return true; - } - private: - /** Compares (operator<) two simplices (i.e. indices) in a - * simplex_order w.r.t. lexicographical order on (value, - * dimension)-tuples. */ - struct cmp_monotone_filtration { - const simplicial_complex &c; - - cmp_monotone_filtration(const simplicial_complex &c) : - c(c){ - } - - bool operator()(index_type i, index_type j) { - const simplex& si = c.simplices[i]; - const simplex& sj = c.simplices[j]; - - if (si.value < sj.value) - return true; - else if (si.value == sj.value) - return si.dim < sj.dim; - else - return false; - } - }; - - public: /** A list of simplices */ std::vector simplices; };