diff options
-rw-r--r-- | src/Simplex_tree/include/gudhi/Simplex_tree.h | 74 | ||||
-rw-r--r-- | src/python/gudhi/simplex_tree.pxd | 2 | ||||
-rw-r--r-- | src/python/gudhi/simplex_tree.pyx | 14 | ||||
-rwxr-xr-x | src/python/test/test_simplex_tree.py | 86 |
4 files changed, 150 insertions, 26 deletions
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 301f7aae..42cf4246 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1467,34 +1467,68 @@ class Simplex_tree { } } - /** \brief Retrieve good values for extended persistence, and separate the diagrams into the ordinary, relative, extended+ and extended- subdiagrams. Need extend_filtration to be called first! - * @param[in] dgm Persistence diagram obtained after calling this->extend_filtration and this->get_persistence. - * @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-. + /** \brief Retrieve good values for extended persistence, and separate the + * diagrams into the ordinary, relative, extended+ and extended- subdiagrams. + * Need extend_filtration to be called first! + * @param[in] dgm Persistence diagram obtained after calling this->extend_filtration + * and this->get_persistence. + * @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-. */ 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; - 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; + std::vector<std::vector<std::pair<int, std::pair<double, double>>>> new_dgm(4); + double x, y; + 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))); + 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. 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. + /** \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. 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. */ void extend_filtration() { // Compute maximum and minimum of filtration values - int maxvert = -std::numeric_limits<double>::infinity(); + int maxvert = -std::numeric_limits<int>::infinity(); std::vector<double> filt; - for (auto sh : this->complex_simplex_range()) {if (this->dimension(sh) == 0){filt.push_back(this->filtration(sh)); maxvert = std::max(*this->simplex_vertex_range(sh).begin(), maxvert);}} + for (auto sh : this->complex_simplex_range()) { + if (this->dimension(sh) == 0){ + filt.push_back(this->filtration(sh)); + maxvert = std::max(*this->simplex_vertex_range(sh).begin(), maxvert); + } + } minval_ = *std::min_element(filt.begin(), filt.end()); maxval_ = *std::max_element(filt.begin(), filt.end()); maxvert += 1; @@ -1502,13 +1536,20 @@ class Simplex_tree { // Compute vectors of integers corresponding to the Simplex handles std::vector<std::vector<int> > splxs; for (auto sh : this->complex_simplex_range()) { - std::vector<int> vr; for (auto vh : this->simplex_vertex_range(sh)){vr.push_back(vh);} + std::vector<int> vr; + for (auto vh : this->simplex_vertex_range(sh)){ + vr.push_back(vh); + } splxs.push_back(vr); } // Add point for coning the simplicial complex int count = this->num_simplices(); - std::vector<int> cone; cone.push_back(maxvert); auto ins = this->insert_simplex(cone, -3); this->assign_key(ins.first, count); count++; + std::vector<int> cone; + cone.push_back(maxvert); + auto ins = this->insert_simplex(cone, -3); + this->assign_key(ins.first, count); + count++; // For each simplex for (auto vr : splxs){ @@ -1531,7 +1572,8 @@ class Simplex_tree { count++; } - this->make_filtration_non_decreasing(); this->initialize_filtration(); + this->make_filtration_non_decreasing(); + this->initialize_filtration(); } diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 4393047f..7aa16926 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -44,7 +44,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": bool prune_above_filtration(double filtration) bool make_filtration_non_decreasing() void extend_filtration() - vector[vector[pair[int, pair[double, double]]]] convert(vector[pair[int, pair[double, double]]]) + vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]]) cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>": diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index cfab14f4..e429e28a 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -387,17 +387,21 @@ cdef class SimplexTree: return self.get_ptr().make_filtration_non_decreasing() def extend_filtration(self): - """ This function extends filtration for computing extended persistence. + """ 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 that after calling this function, the filtration values are actually modified. The function :func:`compute_extended_persistence_subdiagrams()<gudhi.SimplexTree.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. """ return self.get_ptr().extend_filtration() - def convert(self, dgm): - """This function retrieves good values for extended persistence, and separate the diagrams into the ordinary, relative, extended+ and extended- subdiagrams. Need extend_filtration to be called first! + 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 this->extend_filtration and this->get_persistence. + :param dgm: Persistence diagram obtained after calling :func:`extend_filtration()<gudhi.SimplexTree.extend_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-. + + .. note:: + + This function should be called only after calling :func:`extend_filtration()<gudhi.SimplexTree.extend_filtration>` and :func:`persistence()<gudhi.SimplexTree.persistence>`. """ - return self.get_ptr().convert(dgm) + return self.get_ptr().compute_extended_persistence_subdiagrams(dgm) def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 1822c43b..7e3d843e 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -244,7 +244,85 @@ def test_make_filtration_non_decreasing(): assert st.filtration([0, 1, 6]) == 1.0 assert st.filtration([0, 1]) == 1.0 assert st.filtration([0]) == 1.0 - assert st.filtration([1]) == 1.0 - assert st.filtration([3, 4, 5]) == 2.0 - 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 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() + + assert 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))] + ] + |