X-Git-Url: http://git.sthu.org/?a=blobdiff_plain;ds=sidebyside;f=include%2Flibstick%2Fsimplicialcomplex.h;fp=include%2Flibstick%2Fsimplicialcomplex.h;h=1695e25abc6fbe611bb7fa9bfc27e12a618c5f03;hb=44f4198b0d8076203f7247d3183037d5179b11d0;hp=0000000000000000000000000000000000000000;hpb=887a578cf248bc847125debe72549913067b73ba;p=libstick.git diff --git a/include/libstick/simplicialcomplex.h b/include/libstick/simplicialcomplex.h new file mode 100644 index 0000000..1695e25 --- /dev/null +++ b/include/libstick/simplicialcomplex.h @@ -0,0 +1,338 @@ +#ifndef simplicialcomplex_h_nealaezeojeeChuh +#define simplicialcomplex_h_nealaezeojeeChuh + +#include +#include +#include +#include +#include +#include + +#include + +#include "booleanmatrix.h" + + +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. 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 +class simplicial_complex { + + public: + /** The type of this class. */ + 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. */ + 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) { + assert(0 <= dim && dim <= MAXDIM); + + simplex s; + s.dim = dim; + s.value = value; + + if (dim > 0) + 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. */ + 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(); + + return s; + } + + /** Get number of faces. */ + size_t face_count() const { + return face_count_bydim(dim); + } + + /** Get number of faces of a dim-dimensional simplex. */ + static size_t face_count_bydim(int dim) { + assert(-1 <= dim && dim <= MAXDIM); + return dim + 1; + } + }; + + /** An order of the simplices of complex c. An order can be interpreted + * as a permuation of the complex's std::vector of simplices. */ + class simplex_order { + + public: + typedef boolean_colmatrix boundary_matrix; + + /** Create a standard order of the complex c, i.e., the identity permutation. */ + simplex_order(const simplcompltype &c) : + c(c) + { + reset(); + } + + /** Reset order to the identity permutation of the complex's simplices. */ + void reset() { + order.clear(); + for (unsigned i=0; i < c.size(); ++i) + order.push_back(i); + revorder = order; + } + + /** Return number of simplices. */ + size_t size() const { + assert(order.size() == revorder.size()); + return order.size(); + } + + /** Get i-th simplex in the simplex order. */ + 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()); + + 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; + + 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 { + assert(size() == c.size()); + + for (unsigned i=1; i < size(); ++i) + if (get_simplex(i-1).value > get_simplex(i).value) + return false; + + 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()); + + 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. */ + boundary_matrix get_boundary_matrix() const { + 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); + + return mat; + } + + private: + /** Reconstruct 'revorder' by inverting the permutation given by 'order'. */ + void restore_revorder_from_order() { + // Make revorder * order the identity permutation + for (unsigned i=0; i < size(); ++i) + revorder[order[i]] = i; + } + + /** The complex of which we consider a simplex order. */ + const simplcompltype &c; + + /** 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; + + /** 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; + }; + + public: + simplicial_complex() { + // Add the one minus-one dimensional simplex + add_simplex(simplex::create_minusonedim_simplex()); + } + + /** Remove all simplices except the dummy simplex */ + void clear() { + simplices.resize(1); + } + + /** Return number of simplices. */ + size_t size() const { + return simplices.size(); + } + + /** 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); + } + + /** 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(simplex s) { + // Check requirements for faces + for (unsigned i=0; i < s.face_count(); ++i) { + // Faces are already in complex. + assert(s.faces[i] < size()); + // Faces have dimension dim-1 + assert(simplices[s.faces[i]].dim == s.dim-1); + } + + // index_type must be large enough + assert(simplices.size() < std::numeric_limits::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) { + for (unsigned i=0; i < count; ++i) + add_simplex(sarray[i]); + } + + /** 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 { + for (unsigned i=0; i < size(); ++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]]; + if (face.dim != s.dim-1) + return false; + } + } + 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; +}; + +} + + +#endif