Merge pull request #215 from MathieuCarriere/extended_persistence

Extended persistence

4 files changed, 175 insertions, 0 deletions

diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 82f155de..595f22bb 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -57,6 +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 compute_extended_filtration()
+        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 c01cc905..cc3753e1 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -396,6 +396,57 @@ cdef class SimplexTree:
         """
         return self.get_ptr().make_filtration_non_decreasing()
 
+    def extend_filtration(self):
+        """ 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::
+
+            Note that after calling this function, the filtration 
+            values are actually modified within the Simplex_tree. 
+            The function :func:`extended_persistence()<gudhi.SimplexTree.extended_persistence>`
+            retrieves the original values.
+
+        .. 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 possible value (i.e., 4294967295). 
+        """
+        return self.get_ptr().compute_extended_filtration()
+
+    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.
+
+        :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 (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 list of four persistence diagrams in the format described in :func:`persistence()<gudhi.SimplexTree.persistence>`. The first one is Ordinary, the second one is Relative, the third one is Extended+ and the fourth one is Extended-. See https://link.springer.com/article/10.1007/s10208-008-9027-z and/or 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()<gudhi.SimplexTree.extend_filtration>` has been called first!
+
+        .. note::
+
+            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.
+        """
+        cdef vector[pair[int, pair[double, double]]] persistence_result
+        if self.pcohptr != NULL:
+            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):
         """This function returns the persistence of the simplicial complex.
 
diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h
index 4a7062d6..1a18aed6 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<SimplexTreeOptions> {
   using Filtered_simplices = std::vector<Simplex_and_filtration>;
   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,42 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     return cofaces;
   }
 
+  void compute_extended_filtration() {
+    this->efd = this->extend_filtration();
+    this->initialize_filtration();
+    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;
+  }
+
   void create_persistence(Gudhi::Persistent_cohomology_interface<Base>* pcoh) {
     Base::initialize_filtration();
     pcoh = new Gudhi::Persistent_cohomology_interface<Base>(*this);
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index f7848379..70b26e97 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"
@@ -250,6 +251,87 @@ def test_make_filtration_non_decreasing():
     assert st.filtration([3, 4]) == 2.0
     assert st.filtration([4, 5]) == 2.0
 
+def test_extend_filtration():
+
+    # Inserted simplex:
+    #      5   4
+    #      o   o
+    #     / \ /
+    #    o   o
+    #   /2\ /3
+    #  o   o        
+    #  1   0        
+
+    st = SimplexTree()                                                                                                                     
+    st.insert([0,2])
+    st.insert([1,2])
+    st.insert([0,3])
+    st.insert([2,5])
+    st.insert([3,4])
+    st.insert([3,5])                                                                                                         
+    st.assign_filtration([0], 1.)                                                                                                                
+    st.assign_filtration([1], 2.)                                                                                                                
+    st.assign_filtration([2], 3.)                                                                                                                
+    st.assign_filtration([3], 4.)                                                                                                                
+    st.assign_filtration([4], 5.)                                                                                                                
+    st.assign_filtration([5], 6.)                                                                                                                
+
+    assert list(st.get_filtration()) == [                                                                                                                         
+        ([0, 2], 0.0), 
+        ([1, 2], 0.0), 
+        ([0, 3], 0.0), 
+        ([3, 4], 0.0), 
+        ([2, 5], 0.0), 
+        ([3, 5], 0.0), 
+        ([0], 1.0), 
+        ([1], 2.0), 
+        ([2], 3.0), 
+        ([3], 4.0), 
+        ([4], 5.0), 
+        ([5], 6.0)
+    ]
+        
+    st.extend_filtration()
+    
+    assert list(st.get_filtration()) == [                                                                                                                         
+        ([6], -3.0), 
+        ([0], -2.0), 
+        ([1], -1.8), 
+        ([2], -1.6), 
+        ([0, 2], -1.6), 
+        ([1, 2], -1.6), 
+        ([3], -1.4), 
+        ([0, 3], -1.4), 
+        ([4], -1.2), 
+        ([3, 4], -1.2), 
+        ([5], -1.0), 
+        ([2, 5], -1.0), 
+        ([3, 5], -1.0), 
+        ([5, 6], 1.0), 
+        ([4, 6], 1.2), 
+        ([3, 6], 1.4), 
+        ([3, 4, 6], 1.4),
+        ([3, 5, 6], 1.4), 
+        ([2, 6], 1.6), 
+        ([2, 5, 6], 1.6), 
+        ([1, 6], 1.8), 
+        ([1, 2, 6], 1.8), 
+        ([0, 6], 2.0), 
+        ([0, 2, 6], 2.0), 
+        ([0, 3, 6], 2.0)
+    ]
+
+    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.)
+    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():
     st = SimplexTree()