#ifndef image_h_uiHiKahYooroomai
#define image_h_uiHiKahYooroomai
-#include "simplicialcomplex.h"
+#include "simplicialfunction.h"
namespace libstick {
* edges for each of the small squares of the grid graph. Finally the resulting
* triangles are added. For all non-vertex simplices the maximum of its
* vertex-values is used as its value. The first pixel is considered to be at
- * bottom-left. */
-template<class VT, class IT>
-void add_image_to_complex(simplicial_complex<2, IT, VT> &s, const VT *values, size_t width, size_t height) {
+ * bottom-left.
+ *
+ * Requires that f.is_complete() gives true.
+ * */
+template<class IT, class VT>
+void add_image_to_simplicialfunction(simplicial_function<2, IT, VT> &f, const VT *values, size_t width, size_t height) {
+ assert(f.is_complete());
+
// Add the vertices
- const IT v1st = s.size();
+ const IT v1st = f.size();
for (unsigned i=0; i < height; ++i)
for (unsigned j=0; j < width; ++j)
- s.add_simplex(0, NULL, values[i*width + j]);
+ f.add_simplex(0, NULL, values[i*width + j]);
// row i+1 o-------o
// | /|
IT faces[3];
// Adding horizontal edges
- const IT he1st = s.size();
+ const IT he1st = f.size();
for (unsigned i=0; i < height; ++i)
for (unsigned j=0; j < width-1; ++j) {
faces[0] = v1st + i*width + j;
faces[1] = faces[0] + 1;
- s.add_simplex(1, faces);
+ f.add_simplex(1, faces);
}
// Adding vertical edges
- const IT ve1st = s.size();
+ const IT ve1st = f.size();
for (unsigned i=0; i < height-1; ++i)
for (unsigned j=0; j < width; ++j) {
faces[0] = v1st + i*width + j;
faces[1] = faces[0] + width;
- s.add_simplex(1, faces);
+ f.add_simplex(1, faces);
}
// Adding diagonal edges
- const IT de1st = s.size();
+ const IT de1st = f.size();
for (unsigned i=0; i < height-1; ++i)
for (unsigned j=0; j < width-1; ++j) {
faces[0] = v1st + i*width + j;
faces[1] = faces[0] + width + 1;
- s.add_simplex(1, faces);
+ f.add_simplex(1, faces);
}
// Add triangles
faces[0] = he1st + i*(width-1) + j;
faces[1] = ve1st + i*width + j + 1;
faces[2] = de1st + i*(width-1) + j;
- s.add_simplex(2, faces);
+ f.add_simplex(2, faces);
faces[0] = he1st + (i+1)*(width-1) + j;
faces[1] = ve1st + i*width + j;
faces[2] = de1st + i*(width-1) + j;
- s.add_simplex(2, faces);
+ f.add_simplex(2, faces);
}
}
#include <vector>
#include <cassert>
-#include "simplicialcomplex.h"
+#include "simplicialfunction.h"
#include "booleanmatrix.h"
/** Persistence computers various persistent homology information on a
* simplex_order of a simplicial_complex. */
-template<int MAXDIM, class IT, class VT>
+template<int MAXDIM, class IT>
class persistence {
public:
- typedef simplicial_complex<MAXDIM, IT, VT> scomplex;
+ typedef simplicial_complex<MAXDIM, IT> scomplex;
typedef typename scomplex::simplex_order simplex_order;
typedef boolean_colmatrix<IT> boundary_matrix;
typedef boolean_colmatrix<IT> transformation_matrix;
return tm;
}
- /** Print the sequence of births and deaths of cycles (homology classes). */
- void interpret_reduction(std::ostream &os) const {
+ /** Print the sequence of births and deaths of cycles (homology
+ * classes). If a non-NULL simplicial function is given, also the
+ * values of the simplices are printed. */
+ template<class VT>
+ void interpret_reduction(std::ostream &os, simplicial_function<MAXDIM, IT, VT> *f) const {
assert(done_matrices);
for (unsigned c=0; c < bm.width(); ++c) {
os << "simplex";
break;
}
+ os << std::endl;
- os << " \t\tvalue: " << order.get_simplex(c).value << std::endl;
+ if (f != NULL)
+ os << " value: " << f->get_value(order.order_to_complex(c)) << std::endl;
os << " boundary: " << bm.get_column(c) << std::endl;
typename boolean_colmatrix<IT>::column_type col = rbm.get_column(c);
if (col.isempty()) {
os << " \033[1;32mbirth\033[0;m of a cycle: " << tm.get_column(c) << std::endl;
} else {
- os << " \033[1;31mdeath\033[0;m of a cycle: " << rbm.get_column(c) << std::endl;
- os << " boundary of: " << tm.get_column(c) << std::endl;
+ os << " \033[1;31mdeath\033[0;m of a cycle: " << rbm.get_column(c);
+ os << ", boundary of: " << tm.get_column(c) << std::endl;
}
os << std::endl;
}
}
+ /** Print the sequence of births and deaths of cycles (homology classes). */
+ void interpret_reduction(std::ostream &os) const {
+ interpret_reduction<float>(os, NULL);
+ }
+
/** Get the 'dim'-dimensional peristence diagram. */
const diagram& get_persistence_diagram(int dim) const{
assert(-1 <= dim);
#endif
}
- /** Print birth and death information in the persistence diagrams. */
- void interpret_persistence(std::ostream &os) const {
+ /** Print birth and death information in the persistence diagrams. If a
+ * non-NULL simplicial function is given, also the persistence w.r.t.
+ * this function is printed. */
+ template<class VT>
+ void interpret_persistence(std::ostream &os, simplicial_function<MAXDIM, IT, VT> *f) const {
assert(done_diagrams);
for (int d=-1; d <= MAXDIM; ++d) {
os << " ";
if (death > 0) {
- const scomplex& c = get_order().get_complex();
- const VT pers = c.simplices[death].value - c.simplices[birth].value;
- os << birth << "\033[1;31m -> " << death << "\033[0;m" << " \tpers: " << pers;
+ os << birth << "\033[1;31m -> " << death << "\033[0;m";
+
+ if (f != NULL) {
+ const VT pers = f->get_value(death) - f->get_value(birth);
+ os << " \tpers: " << pers;
+ }
} else {
- os << "\033[1;32m" << birth << "\033[0;m" << " \t\t[essential]";
+ os << "\033[1;32m" << birth << "\033[0;m";
+ if (f != NULL)
+ os << " \t\t[essential]";
}
os << std::endl;
}
}
}
+ /** Print birth and death information in the persistence diagrams. */
+ void interpret_persistence(std::ostream &os) const {
+ interpret_persistence<float>(os, NULL);
+ }
+
private:
diagram& _get_persistence_diagram(int dim) {
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;
};
--- /dev/null
+#ifndef simplicialfunction_h_thee5Oonae6eig5k
+#define simplicialfunction_h_thee5Oonae6eig5k
+
+#include <algorithm>
+#include <vector>
+#include <limits>
+
+#include "simplicialcomplex.h"
+
+
+namespace libstick {
+
+/** A simplicial function assigns to each simplex in a
+ * simplicial_complex<MAXDIM, IT> a function value of type VT. A natural way to
+ * obtain a simplex_order is to sort the simplicies of a complex according to
+ * their function values while preserving the filtration property. */
+template<int MAXDIM, class IT, class VT>
+class simplicial_function {
+
+ public:
+ /** The type of the underlying simplicial complex. */
+ typedef simplicial_complex<MAXDIM, IT> simplcompltype;
+ /** Type of indices of simplices. */
+ typedef IT index_type;
+ /** Type of the simplices of the underlying simplicial complex */
+ typedef typename simplcompltype::simplex simplex;
+ /** Type of the order on the underlying simplicial complex */
+ typedef typename simplcompltype::simplex_order simplex_order;
+ /** 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;
+ /** The container that holds the values. */
+ typedef std::vector<value_type> values_type;
+
+ /** This struct adds a function value to a simplex and is designed to
+ * be an aggregate. */
+ struct valuedsimplex {
+ simplex simpl;
+ value_type value;
+
+ /** Create a valuedsimplex. */
+ static valuedsimplex create_valuedsimplex(int dim, index_type* faces, value_type value) {
+ valuedsimplex vs;
+ vs.simpl = simplex::create_simplex(dim, faces);
+ vs.value = value;
+ return vs;
+ }
+ };
+
+
+ /**Create a simplicial function on a given simplicial complex c. Only
+ * define for the (-1)-dim dummy vertex a value, namely
+ * get_minusonedim_value().*/
+ simplicial_function(simplcompltype &c) :
+ c(c) {
+ values.push_back(get_minusonedim_value());
+ }
+
+ /** Get the value of the (-1)-dimensional dummy simplex. */
+ static value_type get_minusonedim_value() {
+ return std::numeric_limits<value_type>::has_infinity
+ ? -std::numeric_limits<value_type>::infinity()
+ : std::numeric_limits<value_type>::min();
+ }
+
+ /** Return number of function values. */
+ size_t size() const {
+ return values.size();
+ }
+
+ /** Returns true if we have a function value for exactly every
+ * simplex. */
+ bool is_complete() const {
+ return (size() == c.size());
+ }
+
+ /** Make the simplicial function complete, i.e. s.t. is_complete()
+ * returns true, by truncating or filling up the list of function
+ * values by the given value. */
+ void make_complete(value_type fill = value_type()) {
+ // We never would like to get rid of the dummy vertex's value
+ assert(c.size() >= 1);
+ values.resize(c.size(), fill);
+ }
+
+ /** Return the underlying complex. */
+ const simplcompltype& get_complex() const {
+ return c;
+ }
+
+ /** Get value of idx-th simplex */
+ const value_type& get_value(index_type idx) const {
+ assert(0 <= idx && idx < size());
+ return values.at(idx);
+ }
+
+ /** Set value of idx-th simplex */
+ void set_value(index_type idx, const value_type& value) {
+ assert(0 <= idx && idx < size());
+ values.at(idx) = value;
+ }
+
+ /** Get all the values. */
+ const values_type& get_values() const {
+ return values;
+ }
+
+ /** Set all the values at once. Requires that is_complete() returns
+ * true afterwards. */
+ void set_values(const values_type& values) {
+ this->values = values;
+ assert(is_complete());
+ }
+
+ /** Add a simplex to the complex and the function. Requires
+ * that is_complete() returns true.
+ *
+ * See also simplicial_complex::add_simplex(). */
+ index_type add_simplex(int dim, index_type* faces, value_type value) {
+ assert(is_complete());
+
+ return add_simplex(valuedsimplex::create_valuedsimplex(dim, faces, value));
+ }
+
+ /** Add a simplex of at least dimension 1 to the complex and the
+ * function. The value of the simplex is set to the maximum value of
+ * its faces. Returns the index of the added simplex. Requires that
+ * is_complete() returns true.
+ *
+ * See also simplicial_complex::add_simplex(). */
+ index_type add_simplex(int dim, index_type* faces) {
+ assert(dim >= 1);
+ assert(is_complete());
+
+ // Get max value of its faces
+ VT value = get_value(faces[0]);
+ for (size_t i=1; i < simplex::face_count_bydim(dim); ++i)
+ value = std::max(value, get_value(faces[i]));
+
+ return add_simplex(dim, faces, value);
+ }
+
+ /** Add a valued simplex. Requires that is_complete() returns true.
+ *
+ * See also simplicial_comple::add_simplex(). */
+ index_type add_simplex(const valuedsimplex &s) {
+ assert(is_complete());
+
+ values.push_back(s.value);
+ return c.add_simplex(s.simpl);
+ }
+
+ /** Add an array of valued simplices. Requires that is_complete()
+ * returns true. */
+ void add_simplices(valuedsimplex* sarray, size_t count) {
+ assert(is_complete());
+
+ for (unsigned i=0; i < count; ++i)
+ add_simplex(sarray[i]);
+ }
+
+ /** 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(size() == c.size());
+ assert(c.is_complex());
+
+ for (unsigned i=0; i < size(); ++i) {
+ const simplex& s = c.get_simplex(i);
+ for (unsigned f=0; f < s.face_count(); ++f)
+ if (get_value(s.faces[f]) > get_value(i))
+ return false;
+ }
+
+ return true;
+ }
+
+ /** Check if the given simplex order is a filtration and monotone
+ * w.r.t. to this simplicial function. Note is_monotone() &&
+ * o.is_filtration() can be true, but is_order_monotonefiltration(o)
+ * may still be false. */
+ bool is_order_monotonefiltration(const simplex_order &o) const {
+ assert(size() == c.size());
+
+ for (unsigned i=1; i < size(); ++i)
+ if (values.at(o.order_to_complex(i-1)) > values.at(o.order_to_complex(i)))
+ return false;
+
+ return o.is_filtration();
+ }
+
+ /** Make simplex order a monotone filtration of the underlying complex. */
+ void make_order_monotonefiltration(simplex_order &so) const {
+ typename simplex_order::order_type order = so.get_order();
+ std::sort(order.begin(), order.end(), cmp_monotone_filtration(*this));
+ so.set_order(order);
+ }
+
+ private:
+ /** Compares (operator<) two valued simplices (i.e. indices) in a
+ * simplicial complex w.r.t. lexicographical order on (value,
+ * dimension)-tuples.
+ * */
+ struct cmp_monotone_filtration {
+ const simplicial_function &f;
+
+ cmp_monotone_filtration(const simplicial_function &f) :
+ f(f) {
+ assert(f.is_complete());
+ }
+
+ bool operator()(index_type i, index_type j) const {
+ assert(f.is_complete());
+
+ if (f.get_value(i) != f.get_value(j))
+ return (f.get_value(i) < f.get_value(j));
+ else
+ return f.get_complex().get_simplex(i).dim < f.get_complex().get_simplex(j).dim;
+ }
+ };
+
+ /** Simplex c[i] get assigned values[i]. */
+ values_type values;
+
+ /** The complex on which we consider a function. */
+ simplcompltype &c;
+};
+
+}
+
+
+#endif
image[i*width + j] = (sin(j/10.0) + 0.5*sin(j/4.0))*cos(i/7.0) + i*j*0.5*1e-3;
// Create complex and add image
- typedef simplicial_complex<2, uint32_t, float> scomplex;
+ typedef simplicial_function<2, uint32_t, float> sfunction;
+ typedef sfunction::simplcompltype scomplex;
scomplex s;
- add_image_to_complex(s, image, width, height);
- assert(s.is_monotone());
+ sfunction f(s);
+ add_image_to_simplicialfunction(f, image, width, height);
+ assert(f.is_monotone());
typedef scomplex::simplex_order sorder;
sorder so(s);
class persistence_TestSuite: public Test::Suite {
private:
- typedef simplicial_complex<3, uint32_t, double> scomplex;
- typedef persistence<3, uint32_t, double> pers;
+ typedef simplicial_function<3, uint32_t, uint32_t> sfunction;
+ typedef sfunction::simplcompltype scomplex;
+ typedef persistence<3, uint32_t> pers;
typedef pers::boundary_matrix bm;
typedef pers::lowestones_matrix lom;
typedef pers::transformation_matrix tm;
// 9
scomplex::simplex ssring[] = {
- // dimension, faces, value...
- /* 1 */ {0, {}, 1},
- /* 2 */ {0, {}, 2},
- /* 3 */ {0, {}, 3},
- /* 4 */ {0, {}, 4},
- /* 5 */ {0, {}, 5},
- /* 6 */ {0, {}, 6},
- /* 7 */ {1, {2, 3}, 6.01},
- /* 8 */ {1, {3, 1}, 6.02},
- /* 9 */ {1, {1, 2}, 6.03},
- /* 10 */ {1, {4, 5}, 6.04},
- /* 11 */ {1, {5, 6}, 6.05},
- /* 12 */ {1, {6, 4}, 6.06},
- /* 13 */ {1, {1, 4}, 6.07},
- /* 14 */ {1, {1, 5}, 6.08},
- /* 15 */ {2, {13, 14, 10}, 6.0801},
- /* 16 */ {1, {3, 6}, 6.09},
- /* 17 */ {1, {2, 5}, 6.10},
- /* 18 */ {1, {2, 6}, 6.11},
- /* 19 */ {2, {9, 14, 17}, 6.1101},
- /* 20 */ {2, {7, 16, 18}, 6.1102},
- /* 21 */ {2, {11, 17, 18}, 6.1103},
- /* 22 */ {1, {3, 4}, 6.12},
- /* 23 */ {2, {12, 16, 22}, 6.1201},
- /* 24 */ {2, {8, 13, 22}, 6.1202}
+ // dimension, faces
+ /* 1 */ {0, {}},
+ /* 2 */ {0, {}},
+ /* 3 */ {0, {}},
+ /* 4 */ {0, {}},
+ /* 5 */ {0, {}},
+ /* 6 */ {0, {}},
+ /* 7 */ {1, {2, 3}},
+ /* 8 */ {1, {3, 1}},
+ /* 9 */ {1, {1, 2}},
+ /* 10 */ {1, {4, 5}},
+ /* 11 */ {1, {5, 6}},
+ /* 12 */ {1, {6, 4}},
+ /* 13 */ {1, {1, 4}},
+ /* 14 */ {1, {1, 5}},
+ /* 15 */ {2, {13, 14, 10}},
+ /* 16 */ {1, {3, 6}},
+ /* 17 */ {1, {2, 5}},
+ /* 18 */ {1, {2, 6}},
+ /* 19 */ {2, {9, 14, 17}},
+ /* 20 */ {2, {7, 16, 18}},
+ /* 21 */ {2, {11, 17, 18}},
+ /* 22 */ {1, {3, 4}},
+ /* 23 */ {2, {12, 16, 22}},
+ /* 24 */ {2, {8, 13, 22}},
};
const size_t cntssring = sizeof(ssring)/sizeof(scomplex::simplex);
ring.add_simplices(ssring, cntssring);
oring.reset();
scomplex::simplex sstorus[] = {
- // dimension, faces, value...
- /* 25 */ {0, {}, 7},
- /* 26 */ {0, {}, 8},
- /* 27 */ {0, {}, 9},
- /* 28 */ {1, {25, 26}, 9.01},
- /* 29 */ {1, {26, 27}, 9.02},
- /* 30 */ {1, {25, 27}, 9.03},
- /* 31 */ {1, {25, 1}, 9.04},
- /* 32 */ {1, {26, 1}, 9.05},
- /* 33 */ {2, {28, 31, 32}, 9.0501},
- /* 34 */ {1, {27, 1}, 9.06},
- /* 35 */ {2, {30, 31, 34}, 9.0601},
- /* 36 */ {1, {27, 3}, 9.07},
- /* 37 */ {2, {36, 34, 8}, 9.0701},
- /* 38 */ {1, {27, 2}, 9.08},
- /* 39 */ {1, {26, 2}, 9.09},
- /* 40 */ {2, {9, 32, 39}, 9.0901},
- /* 41 */ {2, {29, 38, 39}, 9.0902},
- /* 42 */ {2, {7, 38, 36}, 9.0903},
- /* 43 */ {1, {4, 25}, 9.10},
- /* 44 */ {1, {5, 26}, 9.11},
- /* 45 */ {1, {6, 27}, 9.12},
- /* 46 */ {1, {4, 26}, 9.13},
- /* 47 */ {1, {4, 27}, 9.14},
- /* 48 */ {1, {5, 27}, 9.15},
- /* 49 */ {2, {43, 47, 30}, 9.1501},
- /* 50 */ {2, {12, 45, 47}, 9.1502},
- /* 51 */ {2, {29, 44, 48}, 9.1503},
- /* 52 */ {2, {48, 11, 45}, 9.1504},
- /* 53 */ {2, {10, 46, 44}, 9.1505},
- /* 54 */ {2, {43, 46, 28}, 9.1506},
+ // dimension, faces
+ /* 25 */ {0, {}},
+ /* 26 */ {0, {}},
+ /* 27 */ {0, {}},
+ /* 28 */ {1, {25, 26}},
+ /* 29 */ {1, {26, 27}},
+ /* 30 */ {1, {25, 27}},
+ /* 31 */ {1, {25, 1}},
+ /* 32 */ {1, {26, 1}},
+ /* 33 */ {2, {28, 31, 32}},
+ /* 34 */ {1, {27, 1}},
+ /* 35 */ {2, {30, 31, 34}},
+ /* 36 */ {1, {27, 3}},
+ /* 37 */ {2, {36, 34, 8}},
+ /* 38 */ {1, {27, 2}},
+ /* 39 */ {1, {26, 2}},
+ /* 40 */ {2, {9, 32, 39}},
+ /* 41 */ {2, {29, 38, 39}},
+ /* 42 */ {2, {7, 38, 36}},
+ /* 43 */ {1, {4, 25}},
+ /* 44 */ {1, {5, 26}},
+ /* 45 */ {1, {6, 27}},
+ /* 46 */ {1, {4, 26}},
+ /* 47 */ {1, {4, 27}},
+ /* 48 */ {1, {5, 27}},
+ /* 49 */ {2, {43, 47, 30}},
+ /* 50 */ {2, {12, 45, 47}},
+ /* 51 */ {2, {29, 44, 48}},
+ /* 52 */ {2, {48, 11, 45}},
+ /* 53 */ {2, {10, 46, 44}},
+ /* 54 */ {2, {43, 46, 28}},
};
const size_t cntsstorus = sizeof(sstorus)/sizeof(scomplex::simplex);
torus = ring;
otorus.reset();
scomplex::simplex ssball[] = {
- // dimension, faces, value...
- {0, {}, 1},
- {0, {}, 2},
- {0, {}, 3},
- {0, {}, 4},
- {0, {}, 5},
- {1, {1, 2}, 6},
- {1, {2, 3}, 7},
- {1, {3, 4}, 8},
- {1, {4, 1}, 9},
- {1, {1, 5}, 10},
- {1, {2, 5}, 11},
- {1, {3, 5}, 12},
- {1, {4, 5}, 13},
- {2, {6, 10, 11}, 14},
- {2, {7, 11, 12}, 15},
- {2, {8, 12, 13}, 16},
- {2, {9, 13, 10}, 17},
- {1, {1, 3}, 18},
- {2, {18, 6, 7}, 19},
- {2, {18, 8, 9}, 20},
- {2, {18, 10, 12}, 21},
- {3, {21, 14, 15, 19}, 22},
- {3, {21, 16, 17, 20}, 23},
+ // dimension, faces
+ {0, {}},
+ {0, {}},
+ {0, {}},
+ {0, {}},
+ {0, {}},
+ {1, {1, 2}},
+ {1, {2, 3}},
+ {1, {3, 4}},
+ {1, {4, 1}},
+ {1, {1, 5}},
+ {1, {2, 5}},
+ {1, {3, 5}},
+ {1, {4, 5}},
+ {2, {6, 10, 11}},
+ {2, {7, 11, 12}},
+ {2, {8, 12, 13}},
+ {2, {9, 13, 10}},
+ {1, {1, 3}},
+ {2, {18, 6, 7}},
+ {2, {18, 8, 9}},
+ {2, {18, 10, 12}},
+ {3, {21, 14, 15, 19}},
+ {3, {21, 16, 17, 20}},
};
const size_t cntssball = sizeof(ssball)/sizeof(scomplex::simplex);
ball.add_simplices(ssball, cntssball);
lom ballloe(ballloe_coords, ballloe_coords + sizeof(ballloe_coords)/sizeof(uint32_t));
TEST_ASSERT(ballloe == balllo);
- //std::cout << std::endl;
- //ballp.interpret_reduction(std::cout);
- //torusp.compute_diagrams();
- //torusp.interpret_reduction(std::cout);
- //torusp.interpret_persistence(std::cout);
- //std::cout << std::endl;
- //std::cout << std::endl;
+#if 0
+ sfunction ftorus(torus);
+ ftorus.make_complete();
+ for (unsigned i=1; i < torus.size(); ++i)
+ ftorus.set_value(i, i);
+
+ std::cout << std::endl;
+ ballp.interpret_reduction(std::cout);
+ torusp.compute_diagrams();
+ torusp.interpret_reduction(std::cout, &ftorus);
+ torusp.interpret_persistence(std::cout, &ftorus);
+ std::cout << std::endl;
+ std::cout << std::endl;
+#endif
}
void test_betti_numbers() {
class simplicial_complex_TestSuite: public Test::Suite {
private:
- typedef simplicial_complex<3, uint32_t, double> scomplex;
+ typedef simplicial_function<3, uint32_t, double> sfunction;
+ typedef sfunction::simplcompltype scomplex;
typedef scomplex::simplex_order::boundary_matrix bm;
bool setupcalled;
- scomplex c1, c2, c3;
- scomplex::simplex_order o1, o2, o3, o3b;
+ scomplex c;
+ sfunction f1, f2, f3;
+ scomplex::simplex_order o1, o3b;
public:
simplicial_complex_TestSuite() :
setupcalled(false),
- o1(c1),
- o2(c2),
- o3(c3),
- o3b(c3)
+ f1(c),
+ f2(c),
+ f3(c),
+ o1(c),
+ o3b(c)
{
TEST_ADD(simplicial_complex_TestSuite::test_is_complex);
- TEST_ADD(simplicial_complex_TestSuite::test_is_monotoneComplex);
+ TEST_ADD(simplicial_complex_TestSuite::test_is_function_monotone);
TEST_ADD(simplicial_complex_TestSuite::test_is_order_filtration);
TEST_ADD(simplicial_complex_TestSuite::test_is_order_monotone);
TEST_ADD(simplicial_complex_TestSuite::test_boundary_matrix);
return;
setupcalled = true;
- scomplex::simplex ss[] = {
- // dimension, faces, value...
- {0, {}, 1},
- {0, {}, 2},
- {0, {}, 3},
- {0, {}, 4},
- {1, {1, 2}, 5},
- {1, {2, 3}, 6},
- {1, {3, 4}, 7},
- {1, {4, 1}, 8},
- {1, {1, 3}, 9},
- {2, {7, 8, 9}, 10},
- {2, {5, 6, 9}, 11}
+ typedef sfunction::valuedsimplex vsimpl;
+
+ sfunction::valuedsimplex ss[] = {
+ // {dimension, faces}, value
+ {{0, {}}, 1},
+ {{0, {}}, 2},
+ {{0, {}}, 3},
+ {{0, {}}, 4},
+ {{1, {1, 2}}, 5},
+ {{1, {2, 3}}, 6},
+ {{1, {3, 4}}, 7},
+ {{1, {4, 1}}, 8},
+ {{1, {1, 3}}, 9},
+ {{2, {7, 8, 9}}, 10},
+ {{2, {5, 6, 9}}, 11}
};
- const size_t cntss = sizeof(ss)/sizeof(scomplex::simplex);
- c1.add_simplices(ss, cntss);
+ const size_t cntss = sizeof(ss)/sizeof(sfunction::valuedsimplex);
+ f1.add_simplices(ss, cntss);
- // This is o1 This is o2 This is o3(b)
+ // This is f1 This is f2 This is f3
// (value = index) (values) (values)
//
// 1 ----5---- 2 1 ----5---- 2 1 ----5---- 2
// | \| | \| | \|
// 4 ----7---- 3 4 ----7---- 3 4 ----7---- 3
- c2 = c1;
- c2.simplices[6].value = 12;
+ // Copy and change the function values
+ f2.set_values(f1.get_values());
+ f2.set_value(6, 12);
- c3 = c2;
- c3.simplices[11].value = 13;
+ f3.set_values(f2.get_values());
+ f3.set_value(11, 13);
+ // Get two orders, the one is a monotone filtration
o1.reset();
- o2.reset();
- o3.reset();
o3b.reset();
- o3b.make_monotone_filtration();
+ f3.make_order_monotonefiltration(o3b);
}
virtual void tear_down() {
}
void test_is_complex() {
- TEST_ASSERT(c1.is_complex());
- TEST_ASSERT(c2.is_complex());
- TEST_ASSERT(c3.is_complex());
+ TEST_ASSERT(c.is_complex());
}
- void test_is_monotoneComplex() {
- TEST_ASSERT(c1.is_monotone());
- TEST_ASSERT(!c2.is_monotone());
- TEST_ASSERT(c3.is_monotone());
+ void test_is_function_monotone() {
+ TEST_ASSERT(f1.is_monotone());
+ TEST_ASSERT(!f2.is_monotone());
+ TEST_ASSERT(f3.is_monotone());
}
void test_is_order_filtration() {
TEST_ASSERT(o1.is_filtration());
- TEST_ASSERT(o2.is_filtration());
- TEST_ASSERT(o3.is_filtration());
TEST_ASSERT(o3b.is_filtration());
}
void test_is_order_monotone() {
- TEST_ASSERT(o1.is_monotone());
- TEST_ASSERT(!o2.is_monotone());
- TEST_ASSERT(!o3.is_monotone());
- TEST_ASSERT(o3b.is_monotone());
+ TEST_ASSERT(f1.is_order_monotonefiltration(o1));
+ TEST_ASSERT(!f2.is_order_monotonefiltration(o1));
+ TEST_ASSERT(!f3.is_order_monotonefiltration(o1));
+ TEST_ASSERT(f3.is_order_monotonefiltration(o3b));
}
void test_boundary_matrix() {
bm mat1 = o1.get_boundary_matrix();
- bm mat1e(c1.size());
+ bm mat1e(c.size());
uint32_t mat1e_coords[][2] = {
{0, 1}, {0, 2}, {0, 3}, {0, 4}, {1, 5}, {2, 5}, {2, 6}, {3, 6}, {3, 7}, {4, 7}, {1, 8}, {4, 8},
{1, 9}, {3, 9}, {7, 10}, {8, 10}, {9, 10}, {5, 11}, {6, 11}, {9, 11}
bm mat3b = o3b.get_boundary_matrix();
- bm mat3be(c1.size());
+ bm mat3be(c.size());
uint32_t mat3be_coords[][2] = {
{0, 1}, {0, 2}, {0, 3}, {0, 4}, {1, 5}, {2, 5}, {3, 6}, {4, 6}, {1, 7}, {4, 7}, {1, 8}, {3, 8},
{6, 9}, {7, 9}, {8, 9}, {2, 10}, {3, 10}, {5, 11}, {8, 11}, {10, 11}
projects / libstick.git / commitdiff