#include <cstring>
#include <algorithm>
#include <vector>
+#include <limits>
#include <iostream>
-#include <libstick-0.1/booleanmatrix.h>
+#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. */
-template<size_t MAXDIM, class IT=uint32_t, class VT=double>
-class SimplicialComplex {
+ * 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<int MAXDIM, class IT, class VT>
+class simplicial_complex {
public:
/** The type of this class. */
- typedef SimplicialComplex<MAXDIM, IT, VT> simplcompltype;
+ typedef simplicial_complex<MAXDIM, IT, VT> simplcompltype;
/** Type of indices of simplices. */
typedef IT index_type;
/** To every simplex a function value is assigned according to which a
typedef VT value_type;
/** A simplex of the complex. */
- struct Simplex {
+ 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;
+ 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, 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<VT>::has_infinity
+ ? -std::numeric_limits<VT>::infinity()
+ : std::numeric_limits<VT>::min();
return s;
}
}
/** 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 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 SimplexOrder {
+ class simplex_order {
public:
- typedef BooleanColRowMatrix<IT> BoundaryMatrix;
+ typedef boolean_colmatrix<IT> boundary_matrix;
/** Create a standard order of the complex c, i.e., the identity permutation. */
- SimplexOrder(const simplcompltype &c) :
+ simplex_order(const simplcompltype &c) :
c(c)
{
reset();
}
/** Get i-th simplex in the simplex order. */
- const Simplex& getSimplex(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 isFiltration() const {
+ bool is_filtration() const {
assert(size() == c.size());
for (unsigned i=0; i < size(); ++i)
- for (unsigned f=0; f < getSimplex(i).face_count(); ++f)
- if (revorder[getSimplex(i).faces[f]] >= 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 isFiltration() gives true and values of simplices
+ /** Returns true iff is_filtration() gives true and values of simplices
* are monotone w.r.t. this order of simplices. */
- bool isMonotone() const {
+ bool is_monotone() const {
assert(size() == c.size());
for (unsigned i=1; i < size(); ++i)
- if (getSimplex(i-1).value > getSimplex(i).value)
+ if (get_simplex(i-1).value > get_simplex(i).value)
return false;
- return isFiltration();
+ 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 isMonotone() gives true. This
- * requires that the complex's isMonotone() gave true
+ /** Sort simplices such that is_monotone() gives true. This
+ * requires that the complex's is_monotone() gave true
* beforehand.*/
- void makeMonotoneFiltration() {
- assert(c.isMonotone());
+ void make_monotone_filtration() {
+ assert(c.is_monotone());
- sort(order.begin(), order.end(), cmpMonotoneFiltration(c));
- restoreRevorderFromOrder();
+ std::sort(order.begin(), order.end(), cmp_monotone_filtration(c));
+ restore_revorder_from_order();
- assert(c.isMonotone());
- assert(isFiltration());
- assert(isMonotone());
+ assert(c.is_monotone());
+ assert(is_filtration());
+ assert(is_monotone());
}
/** Get the boundary matrix of the complex according to this order. */
- BoundaryMatrix getBoundaryMatrix() const {
- BoundaryMatrix mat(size());
+ boundary_matrix get_boundary_matrix() const {
+ boundary_matrix mat(size());
for (unsigned c=0; c < size(); ++c)
- for(unsigned r=0; r < getSimplex(c).face_count(); ++r)
- mat.set(revorder[getSimplex(c).faces[r]], c, true);
+ 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 restoreRevorderFromOrder() {
+ void restore_revorder_from_order() {
// Make revorder * order the identity permutation
for (unsigned i=0; i < size(); ++i)
revorder[order[i]] = i;
};
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. */
- void addSimplex(size_t dim, index_type* faces, value_type value) {
- addSimplex(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 addSimplex(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.
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) {
+ 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 isComplex() const {
- typename std::vector<Simplex>::const_iterator it = ++simplices.begin();
+ 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;
}
/** 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 isComplex() gives true. */
- bool isMonotone() const {
- assert(isComplex());
+ * the values of its faces. Requires that is_complex() gives true. */
+ 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)
private:
/** Compares (operator<) two simplices (i.e. indices) in a
- * SimplexOrder w.r.t. lexicographical order on (value,
+ * simplex_order w.r.t. lexicographical order on (value,
* dimension)-tuples. */
- struct cmpMonotoneFiltration {
- const SimplicialComplex &c;
+ struct cmp_monotone_filtration {
+ const simplicial_complex &c;
- cmpMonotoneFiltration(const SimplicialComplex &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];
+ const simplex& si = c.simplices[i];
+ const simplex& sj = c.simplices[j];
if (si.value < sj.value)
return true;
public:
/** A list of simplices */
- std::vector<Simplex> simplices;
+ std::vector<simplex> simplices;
};
}
projects / libstick.git / blobdiff