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 +++++++++++++++------------ 1 file changed, 41 insertions(+), 33 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(); } -- cgit v1.2.3