From 7186ae9cc97393496d7f0e0cb281524276f1e47c Mon Sep 17 00:00:00 2001 From: Stefan Huber Date: Thu, 27 Nov 2014 11:47:53 +0100 Subject: [PATCH] Factor out simplicial_function - Factor out simplicial_function from simplical_complex. - simplicial_complex is now independent from the value type VT, and so is the persistence computation. - Simplified test cases for simplicial functions and complexes. --- include/libstick/image.h | 33 ++-- include/libstick/persistence.h | 51 ++++-- include/libstick/simplicialcomplex.h | 230 +++++++++++-------------- include/libstick/simplicialfunction.h | 233 ++++++++++++++++++++++++++ tests/image.h | 8 +- tests/persistence.h | 186 ++++++++++---------- tests/simplicialcomplex.h | 93 +++++----- 7 files changed, 539 insertions(+), 295 deletions(-) create mode 100644 include/libstick/simplicialfunction.h diff --git a/include/libstick/image.h b/include/libstick/image.h index abb1df0..b33801c 100644 --- a/include/libstick/image.h +++ b/include/libstick/image.h @@ -1,7 +1,7 @@ #ifndef image_h_uiHiKahYooroomai #define image_h_uiHiKahYooroomai -#include "simplicialcomplex.h" +#include "simplicialfunction.h" namespace libstick { @@ -14,14 +14,19 @@ 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 -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 +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 // | /| @@ -38,30 +43,30 @@ void add_image_to_complex(simplicial_complex<2, IT, VT> &s, const VT *values, si 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 @@ -70,12 +75,12 @@ void add_image_to_complex(simplicial_complex<2, IT, VT> &s, const VT *values, si 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); } } diff --git a/include/libstick/persistence.h b/include/libstick/persistence.h index b9fc080..a2e163f 100644 --- a/include/libstick/persistence.h +++ b/include/libstick/persistence.h @@ -6,7 +6,7 @@ #include #include -#include "simplicialcomplex.h" +#include "simplicialfunction.h" #include "booleanmatrix.h" @@ -15,11 +15,11 @@ namespace libstick { /** Persistence computers various persistent homology information on a * simplex_order of a simplicial_complex. */ -template +template class persistence { public: - typedef simplicial_complex scomplex; + typedef simplicial_complex scomplex; typedef typename scomplex::simplex_order simplex_order; typedef boolean_colmatrix boundary_matrix; typedef boolean_colmatrix transformation_matrix; @@ -172,8 +172,11 @@ class persistence { 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 + void interpret_reduction(std::ostream &os, simplicial_function *f) const { assert(done_matrices); for (unsigned c=0; c < bm.width(); ++c) { @@ -199,21 +202,28 @@ class persistence { 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::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(os, NULL); + } + /** Get the 'dim'-dimensional peristence diagram. */ const diagram& get_persistence_diagram(int dim) const{ assert(-1 <= dim); @@ -273,8 +283,11 @@ class persistence { #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 + void interpret_persistence(std::ostream &os, simplicial_function *f) const { assert(done_diagrams); for (int d=-1; d <= MAXDIM; ++d) { @@ -287,17 +300,27 @@ class persistence { 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(os, NULL); + } + private: diagram& _get_persistence_diagram(int dim) { 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; }; diff --git a/include/libstick/simplicialfunction.h b/include/libstick/simplicialfunction.h new file mode 100644 index 0000000..6371d66 --- /dev/null +++ b/include/libstick/simplicialfunction.h @@ -0,0 +1,233 @@ +#ifndef simplicialfunction_h_thee5Oonae6eig5k +#define simplicialfunction_h_thee5Oonae6eig5k + +#include +#include +#include + +#include "simplicialcomplex.h" + + +namespace libstick { + +/** A simplicial function assigns to each simplex in a + * simplicial_complex 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 +class simplicial_function { + + public: + /** The type of the underlying simplicial complex. */ + typedef simplicial_complex 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 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::has_infinity + ? -std::numeric_limits::infinity() + : std::numeric_limits::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 diff --git a/tests/image.h b/tests/image.h index a4d7149..9b07eed 100644 --- a/tests/image.h +++ b/tests/image.h @@ -40,10 +40,12 @@ class image_TestSuite: public Test::Suite { 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); diff --git a/tests/persistence.h b/tests/persistence.h index 1028810..186039c 100644 --- a/tests/persistence.h +++ b/tests/persistence.h @@ -12,8 +12,9 @@ using namespace libstick; 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; @@ -64,68 +65,68 @@ class persistence_TestSuite: public Test::Suite { // 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; @@ -133,30 +134,30 @@ class persistence_TestSuite: public Test::Suite { 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); @@ -204,13 +205,20 @@ class persistence_TestSuite: public Test::Suite { 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() { diff --git a/tests/simplicialcomplex.h b/tests/simplicialcomplex.h index da0dfb5..1c8f949 100644 --- a/tests/simplicialcomplex.h +++ b/tests/simplicialcomplex.h @@ -12,23 +12,26 @@ using namespace libstick; 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); @@ -40,24 +43,26 @@ class simplicial_complex_TestSuite: public Test::Suite { 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 @@ -72,51 +77,47 @@ class simplicial_complex_TestSuite: public Test::Suite { // | \| | \| | \| // 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} @@ -126,7 +127,7 @@ class simplicial_complex_TestSuite: public Test::Suite { 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} -- 2.39.5