new fixes

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<Options::store_key, Key_simplex_base_real, Key_simplex_base_dummy>::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<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){
+  std::vector<std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>> extended_persistence_subdiagrams(const std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>& dgm, const Extended_filtration_data& efd){
     std::vector<std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>> 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<Vertex_handle>::min();
-    this->minval_ = std::numeric_limits<Filtration_value>::infinity();
-    this->maxval_ = -std::numeric_limits<Filtration_value>::infinity();
+    Filtration_value minval = std::numeric_limits<Filtration_value>::infinity();
+    Filtration_value maxval = -std::numeric_limits<Filtration_value>::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()<gudhi.SimplexTree.compute_extended_persistence_subdiagrams>`
+            The function :func:`extended_persistence()<gudhi.SimplexTree.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()<gudhi.SimplexTree.extend_filtration>`, :func:`initialize_filtration()<gudhi.SimplexTree.initialize_filtration>`, and :func:`persistence()<gudhi.SimplexTree.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()<gudhi.SimplexTree.extend_filtration>`,
             :func:`initialize_filtration()<gudhi.SimplexTree.initialize_filtration>`, 
-            and :func:`persistence()<gudhi.SimplexTree.persistence>` have been called first!
+            and (optionally) :func:`persistence()<gudhi.SimplexTree.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<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,15 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
     return cofaces;
   }
 
+  void compute_extended_filtration() {
+    this->efd = this->extend_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){
+    return this->extended_persistence_subdiagrams(dgm, this->efd);
+  }
+
   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 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():