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<int MAXDIM, class IT, class VT>
+ * 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<int MAXDIM, class IT>
class simplicial_complex {
public:
/** The type of this class. */
- typedef simplicial_complex<MAXDIM, IT, VT> simplcompltype;
+ typedef simplicial_complex<MAXDIM, IT> 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<VT>::has_infinity
- ? -std::numeric_limits<VT>::infinity()
- : std::numeric_limits<VT>::min();
-
+ simplex s = {-1, {0}};
return s;
}
public:
typedef boolean_colmatrix<IT> boundary_matrix;
+ typedef std::vector<IT> 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();
}
/** 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;
/** 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<index_type> 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<index_type> 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 */
/** 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.
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 {
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<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)
- 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<simplex> simplices;
};
projects / libstick.git / blobdiff