diff options
| -rw-r--r-- | src/Simplex_tree/include/gudhi/Simplex_tree.h | 120 | ||||
| -rw-r--r-- | src/python/gudhi/simplex_tree.pxd | 2 | ||||
| -rw-r--r-- | src/python/gudhi/simplex_tree.pyx | 42 | ||||
| -rwxr-xr-x | src/python/test/test_simplex_tree.py | 83 |
4 files changed, 247 insertions, 0 deletions
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index b7fb9002..1c06e7cb 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -126,6 +126,8 @@ class Simplex_tree { private: typedef typename Dictionary::iterator Dictionary_it; typedef typename Dictionary_it::value_type Dit_value_t; + double minval_; + double maxval_; struct return_first { Vertex_handle operator()(const Dit_value_t& p_sh) const { @@ -1362,6 +1364,7 @@ class Simplex_tree { // Replacing if(f<max) with if(!(f>=max)) would mean that if f is NaN, we replace it with the max of the children. // That seems more useful than keeping NaN. if (!(simplex.second.filtration() >= max_filt_border_value)) { + // Store the filtration modification information modified = true; simplex.second.assign_filtration(max_filt_border_value); @@ -1479,6 +1482,123 @@ class Simplex_tree { } } + /** \brief Retrieve good values for extended persistence, and separate the + * diagrams into the ordinary, relative, extended+ and extended- subdiagrams. + * \post This function should be called only if extend_filtration has been called first! + * \post The coordinates of the persistence diagram points might be a little different than the + * original filtration values due to the internal transformation (scaling to [-2,-1]) that is + * performed on these values during the computation of extended persistence. + * @param[in] dgm Persistence diagram obtained after calling this->extend_filtration, + * this->initialize_filtration, and this->compute_persistent_cohomology. + * @return A vector of four persistence diagrams. The first one is Ordinary, the + * second one is Relative, the third one is Extended+ and the fourth one is Extended-. + * See section 2.2 in https://link.springer.com/article/10.1007/s10208-017-9370-z for a description of these subtypes. + */ + std::vector<std::vector<std::pair<int, std::pair<double, double>>>> compute_extended_persistence_subdiagrams(const std::vector<std::pair<int, std::pair<double, double>>>& dgm){ + std::vector<std::vector<std::pair<int, std::pair<double, double>>>> new_dgm(4); + double x, y; + double minval_ = this->minval_; + double maxval_ = this->maxval_; + for(unsigned int i = 0; i < dgm.size(); i++){ + int h = dgm[i].first; + double px = dgm[i].second.first; + double py = dgm[i].second.second; + if(std::isinf(py)) continue; + else{ + if ((px <= -1) & (py <= -1)){ + x = minval_ + (maxval_-minval_)*(px + 2); + y = minval_ + (maxval_-minval_)*(py + 2); + new_dgm[0].push_back(std::make_pair(h, std::make_pair(x,y))); + } + if ((px >= 1) & (py >= 1)){ + x = minval_ - (maxval_-minval_)*(px - 2); + y = minval_ - (maxval_-minval_)*(py - 2); + new_dgm[1].push_back(std::make_pair(h, std::make_pair(x,y))); + } + if ((px <= -1) & (py >= 1)){ + x = minval_ + (maxval_-minval_)*(px + 2); + y = minval_ - (maxval_-minval_)*(py - 2); + if (x <= y){ + new_dgm[2].push_back(std::make_pair(h, std::make_pair(x,y))); + } + else{ + new_dgm[3].push_back(std::make_pair(h, std::make_pair(x,y))); + } + } + } + } + return new_dgm; + } + + /** \brief Extend filtration for computing extended persistence. + * This function only uses the filtration values at the 0-dimensional simplices, + * and computes the extended persistence diagram induced by the lower-star filtration + * computed with these values. + * \post Note that after calling this function, the filtration + * values are actually modified. The function compute_extended_persistence_subdiagrams + * retrieves the original values and separates the extended persistence diagram points + * w.r.t. their types (Ord, Rel, Ext+, Ext-) and should always be called after + * computing the persistent homology of the extended simplicial complex. + * \post Note that this code creates an extra vertex internally, so you should make sure that + * the Simplex tree does not contain a vertex with the largest Vertex_handle. + */ + void extend_filtration() { + + // Compute maximum and minimum of filtration values + int maxvert = std::numeric_limits<int>::min(); + this->minval_ = std::numeric_limits<double>::max(); + this->maxval_ = std::numeric_limits<double>::min(); + for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh){ + double f = this->filtration(sh); + this->minval_ = std::min(this->minval_, f); + this->maxval_ = std::max(this->maxval_, f); + maxvert = std::max(*this->simplex_vertex_range(sh).begin(), maxvert); + } + + GUDHI_CHECK(maxvert < std::numeric_limits<int>::max(), std::invalid_argument("Simplex_tree contains a vertex with the largest Vertex_handle")); + maxvert += 1; + + Simplex_tree* st_copy = new Simplex_tree(*this); + + // Add point for coning the simplicial complex + int count = this->num_simplices(); + this->insert_simplex({maxvert}, -3); + count++; + + // For each simplex + for (auto sh_copy : st_copy->complex_simplex_range()){ + + // Locate simplex + std::vector<Vertex_handle> vr; + for (auto vh : st_copy->simplex_vertex_range(sh_copy)){ + vr.push_back(vh); + } + auto sh = this->find(vr); + + // Create cone on simplex + vr.push_back(maxvert); + if (this->dimension(sh) == 0){ + // Assign ascending value between -2 and -1 to vertex + double v = this->filtration(sh); + this->assign_filtration(sh, -2 + (v-this->minval_)/(this->maxval_-this->minval_)); + // Assign descending value between 1 and 2 to cone on vertex + this->insert_simplex(vr, 2 - (v-this->minval_)/(this->maxval_-this->minval_)); + } + else{ + // Assign value -3 to simplex and cone on simplex + this->assign_filtration(sh, -3); + this->insert_simplex(vr, -3); + } + count++; + } + + // Deallocate memory + delete st_copy; + + // Automatically assign good values for simplices + this->make_filtration_non_decreasing(); + } + /** \brief Returns a vertex of `sh` that has the same filtration value as `sh` if it exists, and `null_vertex()` otherwise. * * For a lower-star filtration built with `make_filtration_non_decreasing()`, this is a way to invert the process and find out which vertex had its filtration value propagated to `sh`. diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 82f155de..ae32eb82 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -57,6 +57,8 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": void remove_maximal_simplex(vector[int] simplex) bool prune_above_filtration(double filtration) bool make_filtration_non_decreasing() + void extend_filtration() + vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]]) # Iterators over Simplex tree pair[vector[int], double] get_simplex_and_filtration(Simplex_tree_simplex_handle f_simplex) Simplex_tree_simplices_iterator get_simplices_iterator_begin() diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index c01cc905..7af44683 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -396,6 +396,48 @@ cdef class SimplexTree: """ return self.get_ptr().make_filtration_non_decreasing() + def extend_filtration(self): + """ Extend filtration for computing extended persistence. This function only uses the + filtration values at the 0-dimensional simplices, and computes the extended persistence + diagram induced by the lower-star filtration computed with these values. + + .. note:: + + Note that after calling this function, the filtration + values are actually modified within the Simplex_tree. + The function :func:`compute_extended_persistence_subdiagrams()<gudhi.SimplexTree.compute_extended_persistence_subdiagrams>` + retrieves the original values. + + .. note:: + + Note that this code creates an extra vertex internally, so you should make sure that + the Simplex_tree does not contain a vertex with the largest Vertex_handle. + """ + return self.get_ptr().extend_filtration() + + def compute_extended_persistence_subdiagrams(self, dgm): + """This function retrieves good values for extended persistence, and separate the diagrams + into the ordinary, relative, extended+ and extended- subdiagrams. + + :param dgm: Persistence diagram obtained after calling :func:`extend_filtration()<gudhi.SimplexTree.extend_filtration>`, :func:`initialize_filtration()<gudhi.SimplexTree.initialize_filtration>`, and :func:`persistence()<gudhi.SimplexTree.persistence>`. + + :returns: A vector of four persistence diagrams. The first one is Ordinary, the second one is Relative, the third one is Extended+ and the fourth one is Extended-. See section 2.2 in https://link.springer.com/article/10.1007/s10208-017-9370-z for a description of these subtypes. + + .. note:: + + This function should be called only if :func:`extend_filtration()<gudhi.SimplexTree.extend_filtration>`, + :func:`initialize_filtration()<gudhi.SimplexTree.initialize_filtration>`, + and :func:`persistence()<gudhi.SimplexTree.persistence>` have been called first! + + .. note:: + + The coordinates of the persistence diagram points might be a little different than the + original filtration values due to the internal transformation (scaling to [-2,-1]) that is + performed on these values during the computation of extended persistence. + """ + return self.get_ptr().compute_extended_persistence_subdiagrams(dgm) + + def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): """This function returns the persistence of the simplicial complex. diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index f7848379..63eee9a5 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -250,6 +250,89 @@ def test_make_filtration_non_decreasing(): assert st.filtration([3, 4]) == 2.0 assert st.filtration([4, 5]) == 2.0 +def test_extend_filtration(): + + # Inserted simplex: + # 5 4 + # o o + # / \ / + # o o + # /2\ /3 + # o o + # 1 0 + + st = SimplexTree() + st.insert([0,2]) + st.insert([1,2]) + st.insert([0,3]) + st.insert([2,5]) + st.insert([3,4]) + st.insert([3,5]) + st.assign_filtration([0], 1.) + st.assign_filtration([1], 2.) + st.assign_filtration([2], 3.) + st.assign_filtration([3], 4.) + st.assign_filtration([4], 5.) + st.assign_filtration([5], 6.) + + assert list(st.get_filtration()) == [ + ([0, 2], 0.0), + ([1, 2], 0.0), + ([0, 3], 0.0), + ([3, 4], 0.0), + ([2, 5], 0.0), + ([3, 5], 0.0), + ([0], 1.0), + ([1], 2.0), + ([2], 3.0), + ([3], 4.0), + ([4], 5.0), + ([5], 6.0) + ] + + + st.extend_filtration() + st.initialize_filtration() + + assert list(st.get_filtration()) == [ + ([6], -3.0), + ([0], -2.0), + ([1], -1.8), + ([2], -1.6), + ([0, 2], -1.6), + ([1, 2], -1.6), + ([3], -1.4), + ([0, 3], -1.4), + ([4], -1.2), + ([3, 4], -1.2), + ([5], -1.0), + ([2, 5], -1.0), + ([3, 5], -1.0), + ([5, 6], 1.0), + ([4, 6], 1.2), + ([3, 6], 1.4), + ([3, 4, 6], 1.4), + ([3, 5, 6], 1.4), + ([2, 6], 1.6), + ([2, 5, 6], 1.6), + ([1, 6], 1.8), + ([1, 2, 6], 1.8), + ([0, 6], 2.0), + ([0, 2, 6], 2.0), + ([0, 3, 6], 2.0) + ] + + + dgm = st.persistence() + L = st.compute_extended_persistence_subdiagrams(dgm) + assert L == [ + [(0, (1.9999999999999998, 2.9999999999999996))], + [(1, (5.0, 4.0))], + [(0, (1.0, 6.0))], + [(1, (6.0, 1.0))] + ] + + def test_simplices_iterator(): st = SimplexTree() |