From 4c0e6e4144dd3cf6da9600fd4b9bbcac5e664b73 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Sun, 26 Jan 2020 02:54:35 -0500 Subject: added extended persistence function --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 71 +++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) (limited to 'src/Simplex_tree/include/gudhi') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 76608008..4786b244 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -125,6 +125,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 { @@ -1465,6 +1467,75 @@ 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-. + */ + std::vector>>> convert(const std::vector>>& dgm){ + std::vector>>> 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))); + } + } + } + return new_dgm; + } + + /** \brief Extend filtration for computing extended persistence. + */ + void extend_filtration() { + + // Compute maximum and minimum of filtration values + int maxvert = -std::numeric_limits::infinity(); + std::vector 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);}} + minval_ = *std::min_element(filt.begin(), filt.end()); + maxval_ = *std::max_element(filt.begin(), filt.end()); + maxvert += 1; + + // Compute vectors of integers corresponding to the Simplex handles + std::vector > splxs; + for (auto sh : this->complex_simplex_range()) { + std::vector 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 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){ + // Create cone on simplex + auto sh = this->find(vr); 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-minval_)/(maxval_-minval_)); + // Assign descending value between 1 and 2 to cone on vertex + auto ins = this->insert_simplex(vr, 2 - (v-minval_)/(maxval_-minval_)); + this->assign_key(ins.first, count); + } + else{ + // Assign value -3 to simplex and cone on simplex + this->assign_filtration(sh, -3); + auto ins = this->insert_simplex(vr, -3); + this->assign_key(ins.first, count); + } + count++; + } + + this->make_filtration_non_decreasing(); this->initialize_filtration(); + + } + + private: Vertex_handle null_vertex_; /** \brief Total number of simplices in the complex, without the empty simplex.*/ -- cgit v1.2.3 From f2020f6bb3a4d2bbd774aa630151ef1db53ac4f8 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Sun, 2 Feb 2020 15:23:03 -0500 Subject: fixed Marc's comments --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/Simplex_tree/include/gudhi') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 4786b244..301f7aae 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1471,7 +1471,7 @@ class Simplex_tree { * @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>>> convert(const std::vector>>& dgm){ + std::vector>>> compute_extended_persistence_subdiagrams(const std::vector>>& dgm){ std::vector>>> 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; @@ -1487,7 +1487,7 @@ class Simplex_tree { return new_dgm; } - /** \brief Extend filtration for computing extended persistence. + /** \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() { -- cgit v1.2.3 From 360cc2cc31e9e81b99f5c21aa2b4e79b066baabf Mon Sep 17 00:00:00 2001 From: mathieu Date: Tue, 4 Feb 2020 19:44:52 -0500 Subject: fixed Vincent's comments --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 74 ++++++++++++++++++----- src/python/gudhi/simplex_tree.pxd | 2 +- src/python/gudhi/simplex_tree.pyx | 14 +++-- src/python/test/test_simplex_tree.py | 86 +++++++++++++++++++++++++-- 4 files changed, 150 insertions(+), 26 deletions(-) (limited to 'src/Simplex_tree/include/gudhi') 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>>> compute_extended_persistence_subdiagrams(const std::vector>>& dgm){ - std::vector>>> 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>>> 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::infinity(); + int maxvert = -std::numeric_limits::infinity(); std::vector 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 > splxs; for (auto sh : this->complex_simplex_range()) { - std::vector vr; for (auto vh : this->simplex_vertex_range(sh)){vr.push_back(vh);} + std::vector 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 cone; cone.push_back(maxvert); auto ins = this->insert_simplex(cone, -3); this->assign_key(ins.first, count); count++; + std::vector 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>": 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()` 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()` and :func:`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()` and :func:`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))] + ] + -- cgit v1.2.3 From 6e289999fab86bf06cd69c5b7b846c4f26e0a525 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Tue, 17 Mar 2020 00:13:32 -0400 Subject: fixes --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 74 +++++++++++++++------------ src/python/test/test_simplex_tree.py | 12 ++--- 2 files changed, 47 insertions(+), 39 deletions(-) (limited to 'src/Simplex_tree/include/gudhi') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 7be14bce..02f2c7e9 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1354,6 +1354,7 @@ class Simplex_tree { // Replacing 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); @@ -1473,15 +1474,21 @@ 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! + * \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 * 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-. + * See section 2.2 in https://link.springer.com/article/10.1007/s10208-017-9370-z for a description of these subtypes. */ std::vector>>> compute_extended_persistence_subdiagrams(const std::vector>>& dgm){ std::vector>>> 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; @@ -1516,69 +1523,70 @@ class Simplex_tree { /** \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 + * 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::infinity(); - std::vector 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); - } + int maxvert = std::numeric_limits::min(); + this->minval_ = std::numeric_limits::max(); + this->maxval_ = std::numeric_limits::min(); + for (auto sh : this->skeleton_simplex_range(0)) { + 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); } - minval_ = *std::min_element(filt.begin(), filt.end()); - maxval_ = *std::max_element(filt.begin(), filt.end()); + + assert (maxvert < std::numeric_limits::max()); maxvert += 1; - // Compute vectors of integers corresponding to the Simplex handles - std::vector > splxs; - for (auto sh : this->complex_simplex_range()) { - std::vector vr; - for (auto vh : this->simplex_vertex_range(sh)){ - vr.push_back(vh); - } - splxs.push_back(vr); - } + Simplex_tree* st_copy = new Simplex_tree(*this); // Add point for coning the simplicial complex int count = this->num_simplices(); - std::vector cone; - cone.push_back(maxvert); - auto ins = this->insert_simplex(cone, -3); - this->assign_key(ins.first, count); + this->insert_simplex({maxvert}, -3); count++; // For each simplex - for (auto vr : splxs){ + for (auto sh_copy : st_copy->complex_simplex_range()){ + + // Locate simplex + std::vector vr; + for (auto vh : st_copy->simplex_vertex_range(sh_copy)){ + vr.push_back(vh); + } + auto sh = this->find(vr); + // Create cone on simplex - auto sh = this->find(vr); vr.push_back(maxvert); + 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-minval_)/(maxval_-minval_)); + this->assign_filtration(sh, -2 + (v-this->minval_)/(this->maxval_-this->minval_)); // Assign descending value between 1 and 2 to cone on vertex - auto ins = this->insert_simplex(vr, 2 - (v-minval_)/(maxval_-minval_)); - this->assign_key(ins.first, count); + 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); - auto ins = this->insert_simplex(vr, -3); - this->assign_key(ins.first, count); + this->insert_simplex(vr, -3); } count++; } - this->make_filtration_non_decreasing(); - this->initialize_filtration(); + // Deallocate memory + delete st_copy; + // Automatically assign good values for simplices + this->make_filtration_non_decreasing(); } diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index caefeb9c..96ec4707 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -245,6 +245,10 @@ 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(): @@ -271,7 +275,7 @@ def test_extend_filtration(): st.assign_filtration([4], 5.) st.assign_filtration([5], 6.) - assert st.get_filtration() == [ + assert list(st.get_filtration()) == [ ([0, 2], 0.0), ([1, 2], 0.0), ([0, 3], 0.0), @@ -289,7 +293,7 @@ def test_extend_filtration(): st.extend_filtration() - assert st.get_filtration() == [ + assert list(st.get_filtration()) == [ ([6], -3.0), ([0], -2.0), ([1], -1.8), @@ -327,10 +331,6 @@ def test_extend_filtration(): [(1, (6.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_simplices_iterator(): st = SimplexTree() -- cgit v1.2.3 From 58d923b13afb9b18a2d5b028c6575baee691d182 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Tue, 17 Mar 2020 12:14:49 -0400 Subject: update python doc --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 8 +++---- src/python/gudhi/simplex_tree.pyx | 34 +++++++++++++++++++++++---- 2 files changed, 33 insertions(+), 9 deletions(-) (limited to 'src/Simplex_tree/include/gudhi') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 02f2c7e9..f661f687 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1478,8 +1478,8 @@ class Simplex_tree { * \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 - * and this->get_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. @@ -1538,14 +1538,14 @@ class Simplex_tree { int maxvert = std::numeric_limits::min(); this->minval_ = std::numeric_limits::max(); this->maxval_ = std::numeric_limits::min(); - for (auto sh : this->skeleton_simplex_range(0)) { + 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); } - assert (maxvert < std::numeric_limits::max()); + GUDHI_CHECK(maxvert < std::numeric_limits::max(), std::invalid_argument("Simplex_tree contains a vertex with the largest Vertex_handle")); maxvert += 1; Simplex_tree* st_copy = new Simplex_tree(*this); diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 733ecb97..7af44683 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -397,19 +397,43 @@ 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 that after calling this function, the filtration values are actually modified. The function :func:`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. + """ 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()` + 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. + """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()`, :func:`initialize_filtration()`, and :func:`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:: - :param dgm: Persistence diagram obtained after calling :func:`extend_filtration()` and :func:`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-. + This function should be called only if :func:`extend_filtration()`, + :func:`initialize_filtration()`, + and :func:`persistence()` have been called first! .. note:: - This function should be called only after calling :func:`extend_filtration()` and :func:`persistence()`. + 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) -- cgit v1.2.3 From 18a0eb17d9370eca6dde7c0cada0624302ded002 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Tue, 17 Mar 2020 12:31:18 -0400 Subject: implement Marc's suggestions --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) (limited to 'src/Simplex_tree/include/gudhi') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 1c06e7cb..5b36cc1c 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1552,13 +1552,13 @@ class Simplex_tree { 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); + maxvert = std::max(sh->first, maxvert); } GUDHI_CHECK(maxvert < std::numeric_limits::max(), std::invalid_argument("Simplex_tree contains a vertex with the largest Vertex_handle")); maxvert += 1; - Simplex_tree* st_copy = new Simplex_tree(*this); + Simplex_tree st_copy = *this; // Add point for coning the simplicial complex int count = this->num_simplices(); @@ -1566,11 +1566,11 @@ class Simplex_tree { count++; // For each simplex - for (auto sh_copy : st_copy->complex_simplex_range()){ + for (auto sh_copy : st_copy.complex_simplex_range()){ // Locate simplex std::vector vr; - for (auto vh : st_copy->simplex_vertex_range(sh_copy)){ + for (auto vh : st_copy.simplex_vertex_range(sh_copy)){ vr.push_back(vh); } auto sh = this->find(vr); @@ -1592,9 +1592,6 @@ class Simplex_tree { count++; } - // Deallocate memory - delete st_copy; - // Automatically assign good values for simplices this->make_filtration_non_decreasing(); } -- cgit v1.2.3 From a4bf8306d3926428a7d5087d96fbf8033d3bd932 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Tue, 17 Mar 2020 16:17:57 -0400 Subject: fix Marc's comments --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 53 +++++++++++++-------------- 1 file changed, 26 insertions(+), 27 deletions(-) (limited to 'src/Simplex_tree/include/gudhi') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 5b36cc1c..6c837042 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -126,8 +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_; + Filtration_value minval_; + Filtration_value maxval_; struct return_first { Vertex_handle operator()(const Dit_value_t& p_sh) const { @@ -1484,25 +1484,25 @@ 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! + * \pre This function should be called only if this->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. + * @param[in] dgm Persistence diagram obtained after calling this->extend_filtration(), + * this->initialize_filtration(), and Gudhi::persistent_cohomology::Persistent_cohomology::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>>> compute_extended_persistence_subdiagrams(const std::vector>>& dgm){ - std::vector>>> new_dgm(4); - double x, y; - double minval_ = this->minval_; - double maxval_ = this->maxval_; + std::vector>>> compute_extended_persistence_subdiagrams(const std::vector>>& dgm){ + std::vector>>> new_dgm(4); + Filtration_value x, y; + Filtration_value minval_ = this->minval_; + Filtration_value 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; + Filtration_value px = dgm[i].second.first; + Filtration_value py = dgm[i].second.second; if(std::isinf(py)) continue; else{ if ((px <= -1) & (py <= -1)){ @@ -1510,12 +1510,12 @@ class Simplex_tree { 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)){ + else 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)){ + else { x = minval_ + (maxval_-minval_)*(px + 2); y = minval_ - (maxval_-minval_)*(py - 2); if (x <= y){ @@ -1539,37 +1539,36 @@ class Simplex_tree { * 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 + * \pre 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::min(); - this->minval_ = std::numeric_limits::max(); - this->maxval_ = std::numeric_limits::min(); + Vertex_handle maxvert = std::numeric_limits::min(); + this->minval_ = std::numeric_limits::infinity(); + this->maxval_ = -std::numeric_limits::infinity(); for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh){ - double f = this->filtration(sh); + Filtration_value f = this->filtration(sh); this->minval_ = std::min(this->minval_, f); this->maxval_ = std::max(this->maxval_, f); maxvert = std::max(sh->first, maxvert); } - GUDHI_CHECK(maxvert < std::numeric_limits::max(), std::invalid_argument("Simplex_tree contains a vertex with the largest Vertex_handle")); + GUDHI_CHECK(maxvert < std::numeric_limits::max(), std::invalid_argument("Simplex_tree contains a vertex with the largest Vertex_handle")); maxvert += 1; Simplex_tree st_copy = *this; // Add point for coning the simplicial complex - int count = this->num_simplices(); this->insert_simplex({maxvert}, -3); - count++; // For each simplex + std::vector vr; for (auto sh_copy : st_copy.complex_simplex_range()){ // Locate simplex - std::vector vr; + vr.clear(); for (auto vh : st_copy.simplex_vertex_range(sh_copy)){ vr.push_back(vh); } @@ -1578,18 +1577,18 @@ class Simplex_tree { // Create cone on simplex vr.push_back(maxvert); if (this->dimension(sh) == 0){ + Filtration_value v = this->filtration(sh); + Filtration_value scaled_v = (v-this->minval_)/(this->maxval_-this->minval_); // 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_)); + this->assign_filtration(sh, -2 + scaled_v); // Assign descending value between 1 and 2 to cone on vertex - this->insert_simplex(vr, 2 - (v-this->minval_)/(this->maxval_-this->minval_)); + this->insert_simplex(vr, 2 - scaled_v); } else{ // Assign value -3 to simplex and cone on simplex this->assign_filtration(sh, -3); this->insert_simplex(vr, -3); } - count++; } // Automatically assign good values for simplices -- cgit v1.2.3 From 6f445b7e2bdb8481198f8c0f0e076d4fea081d62 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Wed, 18 Mar 2020 12:37:40 -0400 Subject: fix doc --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) (limited to 'src/Simplex_tree/include/gudhi') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 6c837042..697afe26 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1484,12 +1484,12 @@ class Simplex_tree { /** \brief Retrieve good values for extended persistence, and separate the * diagrams into the ordinary, relative, extended+ and extended- subdiagrams. - * \pre This function should be called only if this->extend_filtration() has been called first! + * \pre 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 Gudhi::persistent_cohomology::Persistent_cohomology::compute_persistent_cohomology(). + * @param[in] dgm Persistence diagram obtained after calling `extend_filtration()`, + * `initialize_filtration()`, and `Gudhi::persistent_cohomology::Persistent_cohomology< FilteredComplex, CoefficientField >::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. @@ -1535,7 +1535,7 @@ class Simplex_tree { * 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 + * 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. @@ -1545,7 +1545,7 @@ class Simplex_tree { void extend_filtration() { // Compute maximum and minimum of filtration values - Vertex_handle maxvert = std::numeric_limits::min(); + Vertex_handle maxvert = std::numeric_limits::min(); this->minval_ = std::numeric_limits::infinity(); this->maxval_ = -std::numeric_limits::infinity(); for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh){ -- cgit v1.2.3 From 61691b0081cb868645335c0b1433ddcc0bcbf9e3 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Thu, 19 Mar 2020 13:09:59 -0400 Subject: new fixes --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 45 ++++++++++++++++----------- src/python/gudhi/simplex_tree.pxd | 4 +-- src/python/gudhi/simplex_tree.pyx | 32 ++++++++++++++----- src/python/include/Simplex_tree_interface.h | 13 ++++++++ src/python/test/test_simplex_tree.py | 18 ++++++----- 5 files changed, 77 insertions(+), 35 deletions(-) (limited to 'src/Simplex_tree/include/gudhi') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 697afe26..50b8e582 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -100,6 +100,12 @@ class Simplex_tree { void assign_key(Simplex_key); Simplex_key key() const; }; + struct Extended_filtration_data { + Filtration_value minval; + Filtration_value maxval; + Extended_filtration_data(){} + Extended_filtration_data(Filtration_value vmin, Filtration_value vmax){ minval = vmin; maxval = vmax; } + }; typedef typename std::conditional::type Key_simplex_base; @@ -126,8 +132,6 @@ class Simplex_tree { private: typedef typename Dictionary::iterator Dictionary_it; typedef typename Dictionary_it::value_type Dit_value_t; - Filtration_value minval_; - Filtration_value maxval_; struct return_first { Vertex_handle operator()(const Dit_value_t& p_sh) const { @@ -1490,15 +1494,16 @@ class Simplex_tree { * performed on these values during the computation of extended persistence. * @param[in] dgm Persistence diagram obtained after calling `extend_filtration()`, * `initialize_filtration()`, and `Gudhi::persistent_cohomology::Persistent_cohomology< FilteredComplex, CoefficientField >::compute_persistent_cohomology()`. + * @param[in] efd Structure containing the minimum and maximum values of the original filtration * @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>>> compute_extended_persistence_subdiagrams(const std::vector>>& dgm){ + std::vector>>> extended_persistence_subdiagrams(const std::vector>>& dgm, const Extended_filtration_data& efd){ std::vector>>> new_dgm(4); Filtration_value x, y; - Filtration_value minval_ = this->minval_; - Filtration_value maxval_ = this->maxval_; + Filtration_value minval = efd.minval; + Filtration_value maxval = efd.maxval; for(unsigned int i = 0; i < dgm.size(); i++){ int h = dgm[i].first; Filtration_value px = dgm[i].second.first; @@ -1506,18 +1511,18 @@ class Simplex_tree { if(std::isinf(py)) continue; else{ if ((px <= -1) & (py <= -1)){ - x = minval_ + (maxval_-minval_)*(px + 2); - y = minval_ + (maxval_-minval_)*(py + 2); + 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))); } else if ((px >= 1) & (py >= 1)){ - x = minval_ - (maxval_-minval_)*(px - 2); - y = minval_ - (maxval_-minval_)*(py - 2); + 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))); } else { - x = minval_ + (maxval_-minval_)*(px + 2); - y = minval_ - (maxval_-minval_)*(py - 2); + 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))); } @@ -1535,23 +1540,23 @@ class Simplex_tree { * 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()` + * values are actually modified. The function `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. * \pre 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() { + Extended_filtration_data extend_filtration() { // Compute maximum and minimum of filtration values Vertex_handle maxvert = std::numeric_limits::min(); - this->minval_ = std::numeric_limits::infinity(); - this->maxval_ = -std::numeric_limits::infinity(); + Filtration_value minval = std::numeric_limits::infinity(); + Filtration_value maxval = -std::numeric_limits::infinity(); for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh){ Filtration_value f = this->filtration(sh); - this->minval_ = std::min(this->minval_, f); - this->maxval_ = std::max(this->maxval_, f); + minval = std::min(minval, f); + maxval = std::max(maxval, f); maxvert = std::max(sh->first, maxvert); } @@ -1578,7 +1583,7 @@ class Simplex_tree { vr.push_back(maxvert); if (this->dimension(sh) == 0){ Filtration_value v = this->filtration(sh); - Filtration_value scaled_v = (v-this->minval_)/(this->maxval_-this->minval_); + Filtration_value scaled_v = (v-minval)/(maxval-minval); // Assign ascending value between -2 and -1 to vertex this->assign_filtration(sh, -2 + scaled_v); // Assign descending value between 1 and 2 to cone on vertex @@ -1593,6 +1598,10 @@ class Simplex_tree { // Automatically assign good values for simplices this->make_filtration_non_decreasing(); + + // Return the filtration data + Extended_filtration_data efd(minval, maxval); + return efd; } /** \brief Returns a vertex of `sh` that has the same filtration value as `sh` if it exists, and `null_vertex()` otherwise. diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index ae32eb82..b6284af4 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -57,8 +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]]]) + void compute_extended_filtration() + vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]] dgm) # 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 7af44683..3502000a 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -405,7 +405,7 @@ cdef class SimplexTree: Note that after calling this function, the filtration values are actually modified within the Simplex_tree. - The function :func:`compute_extended_persistence_subdiagrams()` + The function :func:`extended_persistence()` retrieves the original values. .. note:: @@ -413,21 +413,31 @@ cdef class SimplexTree: 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() + return self.get_ptr().compute_extended_filtration() - def compute_extended_persistence_subdiagrams(self, dgm): + def extended_persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): """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()`, :func:`initialize_filtration()`, and :func:`persistence()`. - + :param homology_coeff_field: The homology coefficient field. Must be a + prime number. Default value is 11. + :type homology_coeff_field: int. + :param min_persistence: The minimum persistence value to take into + account (strictly greater than min_persistence). Default value is + 0.0. + Sets min_persistence to -1.0 to see all values. + :type min_persistence: float. + :param persistence_dim_max: If true, the persistent homology for the + maximal dimension in the complex is computed. If false, it is + ignored. Default is false. + :type persistence_dim_max: bool :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()`, :func:`initialize_filtration()`, - and :func:`persistence()` have been called first! + and (optionally) :func:`persistence()` have been called first! .. note:: @@ -435,7 +445,15 @@ cdef class SimplexTree: 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) + cdef vector[pair[int, pair[double, double]]] persistence_result + if self.pcohptr == NULL: + self.pcohptr = new Simplex_tree_persistence_interface(self.get_ptr(), persistence_dim_max) + if self.pcohptr != NULL: + self.pcohptr.get_persistence(homology_coeff_field, min_persistence) + if self.pcohptr != NULL: + pairs = self.pcohptr.persistence_pairs() + persistence_result = [(len(splx1)-1, [self.filtration(splx1), self.filtration(splx2)]) for [splx1, splx2] in pairs] + return self.get_ptr().compute_extended_persistence_subdiagrams(persistence_result) def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 4a7062d6..50ed58d0 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -37,8 +37,12 @@ class Simplex_tree_interface : public Simplex_tree { using Filtered_simplices = std::vector; using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator; using Complex_simplex_iterator = typename Base::Complex_simplex_iterator; + using Extended_filtration_data = typename Base::Extended_filtration_data; public: + + Extended_filtration_data efd; + bool find_simplex(const Simplex& vh) { return (Base::find(vh) != Base::null_simplex()); } @@ -117,6 +121,15 @@ class Simplex_tree_interface : public Simplex_tree { return cofaces; } + void compute_extended_filtration() { + this->efd = this->extend_filtration(); + return; + } + + std::vector>>> compute_extended_persistence_subdiagrams(const std::vector>>& dgm){ + return this->extended_persistence_subdiagrams(dgm, this->efd); + } + void create_persistence(Gudhi::Persistent_cohomology_interface* pcoh) { Base::initialize_filtration(); pcoh = new Gudhi::Persistent_cohomology_interface(*this); diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 63eee9a5..20f6aabf 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -9,6 +9,7 @@ """ from gudhi import SimplexTree +import pytest __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" @@ -322,15 +323,16 @@ def test_extend_filtration(): ([0, 3, 6], 2.0) ] + dgms = st.extended_persistence() - 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))] - ] + assert dgms[0][0][1][0] == pytest.approx(2.) + assert dgms[0][0][1][1] == pytest.approx(3.) + assert dgms[1][0][1][0] == pytest.approx(5.) + assert dgms[1][0][1][1] == pytest.approx(4.) + assert dgms[2][0][1][0] == pytest.approx(1.) + assert dgms[2][0][1][1] == pytest.approx(6.) + assert dgms[3][0][1][0] == pytest.approx(6.) + assert dgms[3][0][1][1] == pytest.approx(1.) def test_simplices_iterator(): -- cgit v1.2.3 From 361abfcfa9ec18c76837f847f8e2e3a060cf7db7 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Thu, 19 Mar 2020 17:02:55 -0400 Subject: added decoding function --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 82 +++++++++++---------------- src/python/gudhi/simplex_tree.pyx | 10 +--- src/python/include/Simplex_tree_interface.h | 27 ++++++++- 3 files changed, 63 insertions(+), 56 deletions(-) (limited to 'src/Simplex_tree/include/gudhi') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 50b8e582..9008c5f2 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -87,6 +87,8 @@ class Simplex_tree { /* \brief Set of nodes sharing a same parent in the simplex tree. */ typedef Simplex_tree_siblings Siblings; + enum Extended_simplex_type {UP, DOWN, EXTRA}; + struct Key_simplex_base_real { Key_simplex_base_real() : key_(-1) {} void assign_key(Simplex_key k) { key_ = k; } @@ -1486,66 +1488,50 @@ class Simplex_tree { } } - /** \brief Retrieve good values for extended persistence, and separate the - * diagrams into the ordinary, relative, extended+ and extended- subdiagrams. + /** \brief Retrieve the original filtration value for a given simplex in the Simplex_tree. Since the + * computation of extended persistence requires modifying the filtration values, this function can be used + * to recover the original values. Moreover, computing extended persistence requires adding new simplices + * in the Simplex_tree. Hence, this function also outputs the type of each simplex. It can be either UP (which means + * that the simplex was present originally, and is thus part of the ascending extended filtration), DOWN (which means + * that the simplex is the cone of an original simplex, and is thus part of the descending extended filtration) or + * EXTRA (which means the simplex is the cone point). Note that if the simplex type is DOWN, the original filtration value + * is set to be the original filtration value of the corresponding (not coned) original simplex. * \pre 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 `extend_filtration()`, - * `initialize_filtration()`, and `Gudhi::persistent_cohomology::Persistent_cohomology< FilteredComplex, CoefficientField >::compute_persistent_cohomology()`. - * @param[in] efd Structure containing the minimum and maximum values of the original filtration - * @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. + * \post The output filtration value is supposed to be the same, but might be a little different, than the + * original filtration value, due to the internal transformation (scaling to [-2,-1]) that is + * performed on the original filtration values during the computation of extended persistence. + * @param[in] f Filtration value of the simplex in the extended (i.e., modified) filtration. + * @param[in] efd Structure containing the minimum and maximum values of the original filtration. This the output of `extend_filtration()`. + * @return A pair containing the original filtration value of the simplex as well as the simplex type. */ - std::vector>>> extended_persistence_subdiagrams(const std::vector>>& dgm, const Extended_filtration_data& efd){ - std::vector>>> new_dgm(4); - Filtration_value x, y; + std::pair decode_extended_filtration(Filtration_value f, const Extended_filtration_data& efd){ + std::pair p; Filtration_value minval = efd.minval; Filtration_value maxval = efd.maxval; - for(unsigned int i = 0; i < dgm.size(); i++){ - int h = dgm[i].first; - Filtration_value px = dgm[i].second.first; - Filtration_value 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))); - } - else 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))); - } - else { - 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 (f >= -2 && f <= -1){ + p.first = minval + (maxval-minval)*(f + 2); p.second = UP; } - return new_dgm; - } + else if (f >= 1 && f <= 2){ + p.first = minval - (maxval-minval)*(f - 2); p.second = DOWN; + } + else{ + p.first = -3; p.second = EXTRA; + } + return p; + }; /** \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 `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. + * values are actually modified. The function `decode_extended_filtration()` + * retrieves the original values and outputs the extended simplex type. * \pre 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. + * the Simplex tree does not contain a vertex with the largest Vertex_handle. + * @return A data structure containing the maximum and minimum values of the original filtration. + * It is meant to be provided as input to `decode_extended_filtration()` in order to retrieve + * the original filtration values for each simplex. */ Extended_filtration_data extend_filtration() { diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 3502000a..2cd81c14 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -415,9 +415,9 @@ cdef class SimplexTree: """ return self.get_ptr().compute_extended_filtration() - def extended_persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): + def extended_persistence(self, homology_coeff_field=11, min_persistence=0): """This function retrieves good values for extended persistence, and separate the diagrams - into the ordinary, relative, extended+ and extended- subdiagrams. + into the Ordinary, Relative, Extended+ and Extended- subdiagrams. :param homology_coeff_field: The homology coefficient field. Must be a prime number. Default value is 11. @@ -427,10 +427,6 @@ cdef class SimplexTree: 0.0. Sets min_persistence to -1.0 to see all values. :type min_persistence: float. - :param persistence_dim_max: If true, the persistent homology for the - maximal dimension in the complex is computed. If false, it is - ignored. Default is false. - :type persistence_dim_max: bool :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:: @@ -447,7 +443,7 @@ cdef class SimplexTree: """ cdef vector[pair[int, pair[double, double]]] persistence_result if self.pcohptr == NULL: - self.pcohptr = new Simplex_tree_persistence_interface(self.get_ptr(), persistence_dim_max) + self.pcohptr = new Simplex_tree_persistence_interface(self.get_ptr(), True) if self.pcohptr != NULL: self.pcohptr.get_persistence(homology_coeff_field, min_persistence) if self.pcohptr != NULL: diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 50ed58d0..a6b1a06e 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -38,6 +38,7 @@ class Simplex_tree_interface : public Simplex_tree { using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator; using Complex_simplex_iterator = typename Base::Complex_simplex_iterator; using Extended_filtration_data = typename Base::Extended_filtration_data; + using Extended_simplex_type = typename Base::Extended_simplex_type; public: @@ -127,7 +128,31 @@ class Simplex_tree_interface : public Simplex_tree { } std::vector>>> compute_extended_persistence_subdiagrams(const std::vector>>& dgm){ - return this->extended_persistence_subdiagrams(dgm, this->efd); + std::vector>>> new_dgm(4); + for (unsigned int i = 0; i < dgm.size(); i++){ + std::pair px = this->decode_extended_filtration(dgm[i].second.first, this->efd); + std::pair py = this->decode_extended_filtration(dgm[i].second.second, this->efd); + std::pair> pd_point = std::make_pair(dgm[i].first, std::make_pair(px.first, py.first)); + //Ordinary + if (px.second == Base::UP && py.second == Base::UP){ + new_dgm[0].push_back(pd_point); + } + // Relative + else if (px.second == Base::DOWN && py.second == Base::DOWN){ + new_dgm[1].push_back(pd_point); + } + else{ + // Extended+ + if (px.first < py.first){ + new_dgm[2].push_back(pd_point); + } + //Extended- + else{ + new_dgm[3].push_back(pd_point); + } + } + } + return new_dgm; } void create_persistence(Gudhi::Persistent_cohomology_interface* pcoh) { -- cgit v1.2.3 From bc223c3cc7cb9e9c0bb3573af720fce9c5380b94 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Mon, 23 Mar 2020 21:22:16 -0400 Subject: new fixes --- src/Simplex_tree/include/gudhi/Simplex_tree.h | 25 +++++++++++++++----- src/python/gudhi/simplex_tree.pxd | 2 +- src/python/gudhi/simplex_tree.pyx | 21 +++++++---------- src/python/include/Simplex_tree_interface.h | 34 ++++++++++++++------------- src/python/test/test_simplex_tree.py | 7 ++---- 5 files changed, 48 insertions(+), 41 deletions(-) (limited to 'src/Simplex_tree/include/gudhi') diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index 9008c5f2..de97d6f2 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -42,6 +42,20 @@ namespace Gudhi { +/** + * \class Extended_simplex_type Simplex_tree.h gudhi/Simplex_tree.h + * \brief Extended simplex type data structure for representing the type of simplices in an extended filtration. + * + * \details The extended simplex type can be either UP (which means + * that the simplex was present originally, and is thus part of the ascending extended filtration), DOWN (which means + * that the simplex is the cone of an original simplex, and is thus part of the descending extended filtration) or + * EXTRA (which means the simplex is the cone point). + * + * Details may be found in section 2.2 in https://link.springer.com/article/10.1007/s10208-017-9370-z. + * + */ +enum class Extended_simplex_type {UP, DOWN, EXTRA}; + struct Simplex_tree_options_full_featured; /** @@ -87,7 +101,7 @@ class Simplex_tree { /* \brief Set of nodes sharing a same parent in the simplex tree. */ typedef Simplex_tree_siblings Siblings; - enum Extended_simplex_type {UP, DOWN, EXTRA}; + struct Key_simplex_base_real { Key_simplex_base_real() : key_(-1) {} @@ -106,7 +120,7 @@ class Simplex_tree { Filtration_value minval; Filtration_value maxval; Extended_filtration_data(){} - Extended_filtration_data(Filtration_value vmin, Filtration_value vmax){ minval = vmin; maxval = vmax; } + Extended_filtration_data(Filtration_value vmin, Filtration_value vmax): minval(vmin), maxval(vmax) {} }; typedef typename std::conditional::type Key_simplex_base; @@ -1370,7 +1384,6 @@ class Simplex_tree { // Replacing 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); @@ -1509,13 +1522,13 @@ class Simplex_tree { Filtration_value minval = efd.minval; Filtration_value maxval = efd.maxval; if (f >= -2 && f <= -1){ - p.first = minval + (maxval-minval)*(f + 2); p.second = UP; + p.first = minval + (maxval-minval)*(f + 2); p.second = Extended_simplex_type::UP; } else if (f >= 1 && f <= 2){ - p.first = minval - (maxval-minval)*(f - 2); p.second = DOWN; + p.first = minval - (maxval-minval)*(f - 2); p.second = Extended_simplex_type::DOWN; } else{ - p.first = -3; p.second = EXTRA; + p.first = std::numeric_limits::quiet_NaN(); p.second = Extended_simplex_type::EXTRA; } return p; }; diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index b6284af4..595f22bb 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -58,7 +58,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": bool prune_above_filtration(double filtration) bool make_filtration_non_decreasing() void compute_extended_filtration() - vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]] dgm) + vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]] dgm, double min_persistence) # 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 5b850462..bcb1578d 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -411,7 +411,7 @@ cdef class SimplexTree: .. 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. + the Simplex_tree does not contain a vertex with the largest possible value (i.e., 4294967295). """ return self.get_ptr().compute_extended_filtration() @@ -422,18 +422,16 @@ cdef class SimplexTree: :param homology_coeff_field: The homology coefficient field. Must be a prime number. Default value is 11. :type homology_coeff_field: int. - :param min_persistence: The minimum persistence value to take into + :param min_persistence: The minimum persistence value (i.e., the absolute value of the difference between the persistence diagram point coordinates) to take into account (strictly greater than min_persistence). Default value is 0.0. Sets min_persistence to -1.0 to see all values. :type min_persistence: float. - :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. + :returns: A list of four persistence diagrams in the format described in :func:`persistence()`. 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()`, - :func:`initialize_filtration()`, - and (optionally) :func:`persistence()` have been called first! + This function should be called only if :func:`extend_filtration()` has been called first! .. note:: @@ -442,14 +440,11 @@ cdef class SimplexTree: performed on these values during the computation of extended persistence. """ cdef vector[pair[int, pair[double, double]]] persistence_result - if self.pcohptr == NULL: - self.pcohptr = new Simplex_tree_persistence_interface(self.get_ptr(), False) - if self.pcohptr != NULL: - self.pcohptr.get_persistence(homology_coeff_field, min_persistence) if self.pcohptr != NULL: - pairs = self.pcohptr.persistence_pairs() - persistence_result = [(len(splx1)-1, [self.filtration(splx1), self.filtration(splx2)]) for [splx1, splx2] in pairs] - return self.get_ptr().compute_extended_persistence_subdiagrams(persistence_result) + del self.pcohptr + self.pcohptr = new Simplex_tree_persistence_interface(self.get_ptr(), False) + persistence_result = self.pcohptr.get_persistence(homology_coeff_field, -1.) + return self.get_ptr().compute_extended_persistence_subdiagrams(persistence_result, min_persistence) def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index a6b1a06e..1a18aed6 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -38,7 +38,6 @@ class Simplex_tree_interface : public Simplex_tree { using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator; using Complex_simplex_iterator = typename Base::Complex_simplex_iterator; using Extended_filtration_data = typename Base::Extended_filtration_data; - using Extended_simplex_type = typename Base::Extended_simplex_type; public: @@ -124,31 +123,34 @@ class Simplex_tree_interface : public Simplex_tree { void compute_extended_filtration() { this->efd = this->extend_filtration(); + this->initialize_filtration(); return; } - std::vector>>> compute_extended_persistence_subdiagrams(const std::vector>>& dgm){ + std::vector>>> compute_extended_persistence_subdiagrams(const std::vector>>& dgm, Filtration_value min_persistence){ std::vector>>> new_dgm(4); for (unsigned int i = 0; i < dgm.size(); i++){ std::pair px = this->decode_extended_filtration(dgm[i].second.first, this->efd); std::pair py = this->decode_extended_filtration(dgm[i].second.second, this->efd); std::pair> pd_point = std::make_pair(dgm[i].first, std::make_pair(px.first, py.first)); - //Ordinary - if (px.second == Base::UP && py.second == Base::UP){ - new_dgm[0].push_back(pd_point); - } - // Relative - else if (px.second == Base::DOWN && py.second == Base::DOWN){ - new_dgm[1].push_back(pd_point); - } - else{ - // Extended+ - if (px.first < py.first){ - new_dgm[2].push_back(pd_point); + if(std::abs(px.first - py.first) > min_persistence){ + //Ordinary + if (px.second == Extended_simplex_type::UP && py.second == Extended_simplex_type::UP){ + new_dgm[0].push_back(pd_point); + } + // Relative + else if (px.second == Extended_simplex_type::DOWN && py.second == Extended_simplex_type::DOWN){ + new_dgm[1].push_back(pd_point); } - //Extended- else{ - new_dgm[3].push_back(pd_point); + // Extended+ + if (px.first < py.first){ + new_dgm[2].push_back(pd_point); + } + //Extended- + else{ + new_dgm[3].push_back(pd_point); + } } } } diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 20f6aabf..70b26e97 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -291,10 +291,8 @@ def test_extend_filtration(): ([5], 6.0) ] - st.extend_filtration() - st.initialize_filtration() - + assert list(st.get_filtration()) == [ ([6], -3.0), ([0], -2.0), @@ -323,7 +321,7 @@ def test_extend_filtration(): ([0, 3, 6], 2.0) ] - dgms = st.extended_persistence() + dgms = st.extended_persistence(min_persistence=-1.) assert dgms[0][0][1][0] == pytest.approx(2.) assert dgms[0][0][1][1] == pytest.approx(3.) @@ -334,7 +332,6 @@ def test_extend_filtration(): assert dgms[3][0][1][0] == pytest.approx(6.) assert dgms[3][0][1][1] == pytest.approx(1.) - def test_simplices_iterator(): st = SimplexTree() -- cgit v1.2.3 From 20ba972d2a7fd14e564ce4adb3921f3f8190fc71 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Wed, 25 Mar 2020 13:00:58 -0400 Subject: update biblio --- biblio/bibliography.bib | 36 +++++++++++++++++++-------- src/Simplex_tree/include/gudhi/Simplex_tree.h | 4 +-- src/python/gudhi/simplex_tree.pyx | 2 +- 3 files changed, 29 insertions(+), 13 deletions(-) (limited to 'src/Simplex_tree/include/gudhi') diff --git a/biblio/bibliography.bib b/biblio/bibliography.bib index 3bbe7960..b017a07e 100644 --- a/biblio/bibliography.bib +++ b/biblio/bibliography.bib @@ -7,11 +7,13 @@ } @article{Carriere17c, - author = {Carri\`ere, Mathieu and Michel, Bertrand and Oudot, Steve}, - title = {{Statistical Analysis and Parameter Selection for Mapper}}, - journal = {CoRR}, - volume = {abs/1706.00204}, - year = {2017} +author = {Carri{\`{e}}re, Mathieu and Michel, Bertrand and Oudot, Steve}, +journal = {Journal of Machine Learning Research}, +pages = {1--39}, +publisher = {JMLR.org}, +title = {{Statistical analysis and parameter selection for Mapper}}, +volume = {19}, +year = {2018} } @inproceedings{Dey13, @@ -23,11 +25,14 @@ } @article{Carriere16, - title={{Structure and Stability of the 1-Dimensional Mapper}}, - author={Carri\`ere, Mathieu and Oudot, Steve}, - journal={CoRR}, - volume= {abs/1511.05823}, - year={2015} +author = {Carri{\`{e}}re, Mathieu and Oudot, Steve}, +journal = {Foundations of Computational Mathematics}, +number = {6}, +pages = {1333--1396}, +publisher = {Springer-Verlag}, +title = {{Structure and stability of the one-dimensional Mapper}}, +volume = {18}, +year = {2017} } @inproceedings{zigzag_reflection, @@ -36,6 +41,17 @@ year = {2014 $\ \ \ \ \ \ \ \ \ \ \ $ \emph{In Preparation}}, } +@article{Cohen-Steiner2009, +author = {Cohen-Steiner, David and Edelsbrunner, Herbert and Harer, John}, +journal = {Foundations of Computational Mathematics}, +number = {1}, +pages = {79--103}, +publisher = {Springer-Verlag}, +title = {{Extending persistence using Poincar{\'{e}} and Lefschetz duality}}, +volume = {9}, +year = {2009} +} + @misc{gudhi_stpcoh, author = {Cl\'ement Maria}, title = "\textsc{Gudhi}, Simplex Tree and Persistent Cohomology Packages", diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index de97d6f2..60720567 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -51,7 +51,7 @@ namespace Gudhi { * that the simplex is the cone of an original simplex, and is thus part of the descending extended filtration) or * EXTRA (which means the simplex is the cone point). * - * Details may be found in section 2.2 in https://link.springer.com/article/10.1007/s10208-017-9370-z. + * Details may be found in \cite Cohen-Steiner2009 and section 2.2 in \cite Carriere16. * */ enum class Extended_simplex_type {UP, DOWN, EXTRA}; @@ -1507,7 +1507,7 @@ class Simplex_tree { * in the Simplex_tree. Hence, this function also outputs the type of each simplex. It can be either UP (which means * that the simplex was present originally, and is thus part of the ascending extended filtration), DOWN (which means * that the simplex is the cone of an original simplex, and is thus part of the descending extended filtration) or - * EXTRA (which means the simplex is the cone point). Note that if the simplex type is DOWN, the original filtration value + * EXTRA (which means the simplex is the cone point). See the definition of Extended_simplex_type. Note that if the simplex type is DOWN, the original filtration value * is set to be the original filtration value of the corresponding (not coned) original simplex. * \pre This function should be called only if `extend_filtration()` has been called first! * \post The output filtration value is supposed to be the same, but might be a little different, than the diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index bcb1578d..6bb22171 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -427,7 +427,7 @@ cdef class SimplexTree: 0.0. Sets min_persistence to -1.0 to see all values. :type min_persistence: float. - :returns: A list of four persistence diagrams in the format described in :func:`persistence()`. 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. + :returns: A list of four persistence diagrams in the format described in :func:`persistence()`. The first one is Ordinary, the second one is Relative, the third one is Extended+ and the fourth one is Extended-. See and https://link.springer.com/article/10.1007/s10208-008-9027-z and section 2.2 in https://link.springer.com/article/10.1007/s10208-017-9370-z for a description of these subtypes. .. note:: -- cgit v1.2.3