From 61691b0081cb868645335c0b1433ddcc0bcbf9e3 Mon Sep 17 00:00:00 2001 From: MathieuCarriere Date: Thu, 19 Mar 2020 13:09:59 -0400 Subject: new fixes --- src/python/include/Simplex_tree_interface.h | 13 +++++++++++++ 1 file changed, 13 insertions(+) (limited to 'src/python/include/Simplex_tree_interface.h') 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); -- 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/python/include/Simplex_tree_interface.h') 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/python/include/Simplex_tree_interface.h') 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