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(-) 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