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 +++++++++++++++++++++++++++ src/python/gudhi/simplex_tree.pxd | 2 + src/python/gudhi/simplex_tree.pyx | 14 ++++++ 3 files changed, 87 insertions(+) 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.*/ diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 1066d44b..39f2a45f 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -43,6 +43,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]]]] convert(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 b18627c4..cfab14f4 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -386,6 +386,20 @@ cdef class SimplexTree: """ return self.get_ptr().make_filtration_non_decreasing() + def extend_filtration(self): + """ This function extends filtration for computing extended persistence. + """ + 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! + + :param dgm: Persistence diagram obtained after calling this->extend_filtration and this->get_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-. + """ + return self.get_ptr().convert(dgm) + + def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): """This function returns the persistence of the simplicial complex. -- cgit v1.2.3