extended_persistence uses directly get_persistent_pairs

5 files changed, 47 insertions, 34 deletions

diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 5c98fb4a..a8ed6d50 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -63,7 +63,6 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         bool prune_above_filtration(double filtration) nogil
         bool make_filtration_non_decreasing() nogil
         void compute_extended_filtration() nogil
-        vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]] dgm, double min_persistence) nogil
         Simplex_tree_interface_full_featured* collapse_edges(int nb_collapse_iteration) nogil except +
         void reset_filtration(double filtration, int dimension) nogil
         bint operator==(Simplex_tree_interface_full_featured) nogil
@@ -78,7 +77,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         pair[Simplex_tree_boundary_iterator, Simplex_tree_boundary_iterator] get_boundary_iterators(vector[int] simplex) nogil except +
 
 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
-    cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>":
+    cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>>":
         Simplex_tree_persistence_interface(Simplex_tree_interface_full_featured * st, bool persistence_dim_max) nogil
         void compute_persistence(int homology_coeff_field, double min_persistence) nogil except +
         vector[pair[int, pair[double, double]]] get_persistence() nogil
@@ -89,3 +88,4 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
         vector[pair[vector[int], vector[int]]] persistence_pairs() nogil
         pair[vector[vector[int]], vector[vector[int]]] lower_star_generators() nogil
         pair[vector[vector[int]], vector[vector[int]]] flag_generators() nogil
+        vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(double min_persistence) nogil
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index b8fabf78..1bbf1539 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -471,8 +471,7 @@ cdef class SimplexTree:
             del self.pcohptr
         self.pcohptr = new Simplex_tree_persistence_interface(self.get_ptr(), False)
         self.pcohptr.compute_persistence(homology_coeff_field, -1.)
-        persistence_result = self.pcohptr.get_persistence()
-        return self.get_ptr().compute_extended_persistence_subdiagrams(persistence_result, min_persistence)
+        return self.pcohptr.compute_extended_persistence_subdiagrams(min_persistence)
 
 
     def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False):
diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h
index e5a3dfba..945378a0 100644
--- a/src/python/include/Persistent_cohomology_interface.h
+++ b/src/python/include/Persistent_cohomology_interface.h
@@ -12,6 +12,8 @@
 #define INCLUDE_PERSISTENT_COHOMOLOGY_INTERFACE_H_
 
 #include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Simplex_tree.h>  // for Extended_simplex_type
+
 
 #include <cstdlib>
 #include <vector>
@@ -223,6 +225,44 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol
     return out;
   }
 
+  using Filtration_value = typename FilteredComplex::Filtration_value;
+  using Birth_death = std::pair<Filtration_value, Filtration_value>;
+  using Persistence_subdiagrams = std::vector<std::vector<std::pair<int, Birth_death>>>;
+
+  Persistence_subdiagrams compute_extended_persistence_subdiagrams(Filtration_value min_persistence){
+    Persistence_subdiagrams pers_subs(4);
+    auto const& persistent_pairs = Base::get_persistent_pairs();
+    for (auto pair : persistent_pairs) {
+      std::pair<Filtration_value, Extended_simplex_type> px = stptr_->decode_extended_filtration(stptr_->filtration(get<0>(pair)),
+                                                                                                        stptr_->efd);
+      std::pair<Filtration_value, Extended_simplex_type> py = stptr_->decode_extended_filtration(stptr_->filtration(get<1>(pair)),
+                                                                                                        stptr_->efd);
+      std::pair<int, Birth_death> pd_point = std::make_pair(stptr_->dimension(get<0>(pair)),
+                                                            std::make_pair(px.first, py.first));
+      if(std::abs(px.first - py.first) > min_persistence){
+        //Ordinary
+        if (px.second == Extended_simplex_type::UP && py.second == Extended_simplex_type::UP){
+          pers_subs[0].push_back(pd_point);
+        }
+        // Relative
+        else if (px.second == Extended_simplex_type::DOWN && py.second == Extended_simplex_type::DOWN){
+          pers_subs[1].push_back(pd_point);
+        }
+        else{
+          // Extended+
+          if (px.first < py.first){
+            pers_subs[2].push_back(pd_point);
+          }
+          //Extended-
+          else{
+            pers_subs[3].push_back(pd_point);
+          }
+        }
+      }
+    }
+    return pers_subs;
+  }
+
  private:
   // A copy
   FilteredComplex* stptr_;
diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h
index 629f6083..dc9d01d7 100644
--- a/src/python/include/Simplex_tree_interface.h
+++ b/src/python/include/Simplex_tree_interface.h
@@ -133,36 +133,6 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     return;
   }
 
-  std::vector<std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>> compute_extended_persistence_subdiagrams(const std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>& dgm, Filtration_value min_persistence){
-    std::vector<std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>> new_dgm(4);
-    for (unsigned int i = 0; i < dgm.size(); i++){
-      std::pair<Filtration_value, Extended_simplex_type> px = this->decode_extended_filtration(dgm[i].second.first, this->efd);
-      std::pair<Filtration_value, Extended_simplex_type> py = this->decode_extended_filtration(dgm[i].second.second, this->efd);
-      std::pair<int, std::pair<Filtration_value, Filtration_value>> pd_point = std::make_pair(dgm[i].first, std::make_pair(px.first, py.first));
-      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);
-        }
-        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;
-  }
-
   Simplex_tree_interface* collapse_edges(int nb_collapse_iteration) {
 #ifdef GUDHI_USE_EIGEN3
     using Filtered_edge = std::tuple<Vertex_handle, Vertex_handle, Filtration_value>;
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index a8180ce8..23458eb2 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -320,6 +320,10 @@ def test_extend_filtration():
     ]
 
     dgms = st.extended_persistence(min_persistence=-1.)
+    assert len(dgms) == 4
+    # Sort by (death-birth) descending - we are only interested in those with the longest life span
+    for idx in range(4):
+        dgms[idx] = sorted(dgms[idx], key=lambda x:(-abs(x[1][0]-x[1][1])))
 
     assert dgms[0][0][1][0] == pytest.approx(2.)
     assert dgms[0][0][1][1] == pytest.approx(3.)