fix Marc's comments

2 files changed, 55 insertions, 61 deletions

diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx
index bd432834..8cf43539 100644
--- a/src/python/gudhi/cubical_complex.pyx
+++ b/src/python/gudhi/cubical_complex.pyx
@@ -31,7 +31,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
     cdef cppclass Cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Cubical_complex::Cubical_complex_interface<>>":
         Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max)
         vector[pair[int, pair[double, double]]] get_persistence(int homology_coeff_field, double min_persistence)
-        vector[pair[int, pair[pair[double, int], pair[double, int]]]] get_cofaces_of_cubical_persistence_pairs(int homology_coeff_field, double min_persistence)
+        vector[vector[int]] cofaces_of_cubical_persistence_pairs()
         vector[int] betti_numbers()
         vector[int] persistent_betti_numbers(double from_value, double to_value)
         vector[pair[double,double]] intervals_in_dimension(int dimension)
@@ -146,31 +146,32 @@ cdef class CubicalComplex:
             persistence_result = self.pcohptr.get_persistence(homology_coeff_field, min_persistence)
         return persistence_result
 
-    def cofaces_of_cubical_persistence_pairs(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False):
-        """This function returns the persistence of the simplicial complex.
-
-        :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: The persistence of the simplicial complex, together with the cofaces of the corresponding generators, i.e., the top-dimensional cells/cofaces of the positive and negative simplices.
-        :rtype:  list of pairs(dimension, pair(index of positive top-dimensional cell, index of negative top-dimensional cell))
+    def cofaces_of_persistence_pairs(self):
+        """A persistence interval is described by a pair of cells, one that creates the 
+           feature and one that kills it. The filtration values of those 2 cells give coordinates 
+           for a point in a persistence diagram, or a bar in a barcode. Structurally, in the 
+           cubical complexes provided here, the filtration value of any cell is the minimum of the 
+           filtration values of the maximal cells that contain it. Connecting persistence diagram 
+           coordinates to the corresponding value in the input (i.e. the filtration values of 
+           the top-dimensional cells) is useful for differentiation purposes.
+
+           This function returns a list of pairs of top-dimensional cells corresponding to 
+           the persistence birth and death cells of the filtration. The cells are represented by 
+           their indices in the input list of top-dimensional cells (and not their indices in the 
+           internal datastructure that includes non-maximal cells). Note that when two adjacent 
+           top-dimensional cells have the same filtration value, we arbitrarily return one of the two
+           when calling the function on one of their common faces.
+
+        :returns: The top-dimensional cells/cofaces of the positive and negative cells.
+        :rtype:  list of pairs(index of positive top-dimensional cell, index of negative top-dimensional cell)
         """
+        cdef vector[vector[int]] persistence_result
         if self.pcohptr != NULL:
-            del self.pcohptr
-        self.pcohptr = new Cubical_complex_persistence_interface(self.thisptr, True)
-        cdef vector[pair[int, pair[pair[double, int], pair[double, int]]]] persistence_result
-        if self.pcohptr != NULL:
-            persistence_result = self.pcohptr.get_cofaces_of_cubical_persistence_pairs(homology_coeff_field, min_persistence)
-        return persistence_result
+            persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs()
+        else:
+            print("cofaces_of_persistence_pairs function requires persistence function"
+                  " to be launched first.")
+        return np.array(persistence_result)
 
     def betti_numbers(self):
         """This function returns the Betti numbers of the complex.
diff --git a/src/python/include/Persistent_cohomology_interface.h b/src/python/include/Persistent_cohomology_interface.h
index 1a1e716e..e5accf50 100644
--- a/src/python/include/Persistent_cohomology_interface.h
+++ b/src/python/include/Persistent_cohomology_interface.h
@@ -63,7 +63,6 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol
     auto persistent_pairs = persistent_cohomology::Persistent_cohomology<FilteredComplex,
       persistent_cohomology::Field_Zp>::get_persistent_pairs();
     std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp);
-
     std::vector<std::pair<int, std::pair<double, double>>> persistence;
     for (auto pair : persistent_pairs) {
       persistence.push_back(std::make_pair(stptr_->dimension(get<0>(pair)),
@@ -73,58 +72,52 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol
     return persistence;
   }
 
-  void top_dimensional_cofaces(std::vector<int> & cofaces, int splx){
-    if (stptr_->dimension(stptr_->simplex(splx)) == stptr_->dimension()){cofaces.push_back(stptr_->simplex(splx));}
-    else{  for (auto v : stptr_->coboundary_simplex_range(stptr_->simplex(splx))){top_dimensional_cofaces(cofaces, stptr_->key(v));}  }
+  int top_dimensional_coface(int splx){
+    if (stptr_->dimension(splx) == stptr_->dimension()){return splx;}
+    else{  
+      for (auto v : stptr_->coboundary_simplex_range(splx)){
+        if(stptr_->filtration(v) == stptr_->filtration(splx)){
+          return top_dimensional_coface(v);
+        }
+      }  
+    }
   }
 
-  std::vector<std::pair<int, std::pair<std::pair<double, int>, std::pair<double, int>>>> get_cofaces_of_cubical_persistence_pairs(int homology_coeff_field,
-                                                                                                                                  double min_persistence) {
+  std::vector<std::vector<int>> cofaces_of_cubical_persistence_pairs() {
 
     // Warning: this function is meant to be used with CubicalComplex only!!
 
+    auto pairs = persistent_cohomology::Persistent_cohomology<FilteredComplex,
+      persistent_cohomology::Field_Zp>::get_persistent_pairs();
+
     // Gather all top-dimensional cells and store their simplex handles
-    std::vector<int> max_splx; for (auto splx : stptr_->filtration_simplex_range()){ if (stptr_->dimension(splx) == stptr_->dimension())  max_splx.push_back(splx); }
+    std::vector<int> max_splx; for (auto splx : stptr_->top_dimensional_cells_range()){
+      max_splx.push_back(splx); 
+    }
     // Sort these simplex handles and compute the ordering function
     // This function allows to go directly from the simplex handle to the position of the corresponding top-dimensional cell in the input data 
-    std::map<int, int> order; std::sort(max_splx.begin(), max_splx.end()); for (unsigned int i = 0; i < max_splx.size(); i++)  order.insert(std::make_pair(max_splx[i], i));
-
-    persistent_cohomology::Persistent_cohomology<FilteredComplex,
-      persistent_cohomology::Field_Zp>::init_coefficients(homology_coeff_field);
-    persistent_cohomology::Persistent_cohomology<FilteredComplex,
-      persistent_cohomology::Field_Zp>::compute_persistent_cohomology(min_persistence);
+    std::map<int, int> order; std::sort(max_splx.begin(), max_splx.end()); 
+    for (unsigned int i = 0; i < max_splx.size(); i++)  order.insert(std::make_pair(max_splx[i], i));
 
-    // Custom sort and output persistence
-    cmp_intervals_by_dim_then_length cmp(stptr_);
-    auto persistent_pairs = persistent_cohomology::Persistent_cohomology<FilteredComplex,
-      persistent_cohomology::Field_Zp>::get_persistent_pairs();
-    std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp);
-
-    std::vector<std::pair<int, std::pair<std::pair<double, int>, std::pair<double, int>>>> persistence;
-    for (auto pair : persistent_pairs) {
-
-      double f0 = stptr_->filtration(get<0>(pair));
-      // Recursively get the top-dimensional cells / cofaces associated to the persistence generator
-      std::vector<int> faces0; top_dimensional_cofaces(faces0, stptr_->key(get<0>(pair)));
-      // Find the top-dimensional cell / coface with the same filtration value
-      int cf; for (unsigned int i = 0; i < faces0.size(); i++){if (stptr_->filtration(faces0[i]) == f0){cf = i; break;}}
+    std::vector<std::vector<int>> persistence_pairs;
+    for (auto pair : pairs) {
+      int h = stptr_->dimension(get<0>(pair));
+      // Recursively get the top-dimensional cell / coface associated to the persistence generator
+      int face0 = top_dimensional_coface(get<0>(pair));
       // Retrieve the index of the corresponding top-dimensional cell in the input data
-      int splx0 = order[faces0[cf]];
+      int splx0 = order[face0];
 
       int splx1 = -1;
       if (isfinite(stptr_->filtration(get<1>(pair)))){
-      double f1 = stptr_->filtration(get<1>(pair));
-      // Recursively get the top-dimensional cells / cofaces associated to the persistence generator
-      std::vector<int> faces1; top_dimensional_cofaces(faces1, stptr_->key(get<1>(pair)));
-      // Find the top-dimensional cell / coface with the same filtration value
-      int cf; for (unsigned int i = 0; i < faces1.size(); i++){if (stptr_->filtration(faces1[i]) == f1){cf = i; break;}}
+      // Recursively get the top-dimensional cell / coface associated to the persistence generator
+      int face1 = top_dimensional_coface(get<1>(pair));
       // Retrieve the index of the corresponding top-dimensional cell in the input data
-      splx1 = order[faces1[cf]];
+      splx1 = order[face1];
       }
-
-      persistence.push_back(std::make_pair(stptr_->dimension(get<0>(pair)), std::make_pair(std::make_pair(stptr_->filtration(get<0>(pair)), splx0), std::make_pair(stptr_->filtration(get<1>(pair)), splx1))));
+      std::vector<int> vect{ h, splx0, splx1};
+      persistence_pairs.push_back(vect);
     }
-    return persistence;
+    return persistence_pairs;
   }
 
   std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs() {