update output

5 files changed, 69 insertions, 15 deletions

diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
index 6441c129..248ebdb6 100644
--- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
+++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex_base.h
@@ -110,8 +110,9 @@ class Bitmap_cubical_complex_base {
   virtual inline std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const;
 
   /**
-   * This function computes the index of one of the top-dimensional cubes (chosen arbitrarily) associated
-   * to a given simplex handle. Note that the input parameter is not necessarily a cube, it might also
+   * This function finds a top-dimensional cell that is incident to the input cell and has 
+   * the same filtration value. In case several cells are suitable, an arbitrary one is 
+   * returned. Note that the input parameter is not necessarily a cube, it might also
    * be an edge or vertex of a cube. On the other hand, the output is always indicating the position of
    * a cube in the data structure. 
    **/
diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx
index 69d0f0b6..884b0664 100644
--- a/src/python/gudhi/cubical_complex.pyx
+++ b/src/python/gudhi/cubical_complex.pyx
@@ -187,18 +187,41 @@ cdef class CubicalComplex:
         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, together with the corresponding homological dimension.
-        :rtype:  numpy array of integers of shape [number_of_persistence_points, 3], the integers of eah row being: (homological dimension, 
-            index of positive top-dimensional cell, index of negative top-dimensional cell). If the homological feature is essential, i.e., if 
-            the death time is +infinity, then the index of the corresponding negative top-dimensional cell is -1.
+        :returns: The top-dimensional cells/cofaces of the positive and negative cells, 
+            together with the corresponding homological dimension, in two lists of numpy arrays of integers.
+            The first list contains the regular persistence pairs, grouped by dimension. 
+            It contains numpy arrays of shape [number_of_persistence_points, 2].
+            The indices of the arrays in the list correspond to the homological dimensions, and the 
+            integers of each row in each array correspond to: (index of positive top-dimensional cell, 
+            index of negative top-dimensional cell). 
+            The second list contains the essential features, grouped by dimension.             
+            It contains numpy arrays of shape [number_of_persistence_points, 1].
+            The indices of the arrays in the list correspond to the homological dimensions, and the 
+            integers of each row in each array correspond to: (index of positive top-dimensional cell).
         """
         cdef vector[vector[int]] persistence_result
         if self.pcohptr != NULL:
+            output = [[],[]]
             persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs()
+            pr = np.array(persistence_result)
+
+            ess_ind = np.argwhere(pr[:,2] == -1)[:,0]
+            ess = pr[ess_ind]
+            max_h = max(ess[:,0])+1
+            for h in range(max_h):
+                hidxs = np.argwhere(ess[:,0] == h)[:,0]
+                output[1].append(ess[hidxs][:,1])
+
+            reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind)
+            reg = pr[reg_ind]
+            max_h = max(reg[:,0])+1
+            for h in range(max_h):
+                hidxs = np.argwhere(reg[:,0] == h)[:,0]
+                output[0].append(reg[hidxs][:,1:])
         else:
             print("cofaces_of_persistence_pairs function requires persistence function"
                   " to be launched first.")
-        return np.array(persistence_result)
+        return output
 
     def betti_numbers(self):
         """This function returns the Betti numbers of the complex.
diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx
index 78565cf8..3cf2ff01 100644
--- a/src/python/gudhi/periodic_cubical_complex.pyx
+++ b/src/python/gudhi/periodic_cubical_complex.pyx
@@ -192,18 +192,41 @@ cdef class PeriodicCubicalComplex:
         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, together with the corresponding homological dimension.
-        :rtype:  numpy array of integers of shape [number_of_persistence_points, 3], the integers of eah row being: (homological dimension, 
-            index of positive top-dimensional cell, index of negative top-dimensional cell). If the homological feature is essential, i.e., if 
-            the death time is +infinity, then the index of the corresponding negative top-dimensional cell is -1.
+        :returns: The top-dimensional cells/cofaces of the positive and negative cells, 
+            together with the corresponding homological dimension, in two lists of numpy arrays of integers.
+            The first list contains the regular persistence pairs, grouped by dimension. 
+            It contains numpy arrays of shape [number_of_persistence_points, 2].
+            The indices of the arrays in the list correspond to the homological dimensions, and the 
+            integers of each row in each array correspond to: (index of positive top-dimensional cell, 
+            index of negative top-dimensional cell). 
+            The second list contains the essential features, grouped by dimension.             
+            It contains numpy arrays of shape [number_of_persistence_points, 1].
+            The indices of the arrays in the list correspond to the homological dimensions, and the 
+            integers of each row in each array correspond to: (index of positive top-dimensional cell).
         """
         cdef vector[vector[int]] persistence_result
         if self.pcohptr != NULL:
+            output = [[],[]]
             persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs()
+            pr = np.array(persistence_result)
+
+            ess_ind = np.argwhere(pr[:,2] == -1)[:,0]
+            ess = pr[ess_ind]
+            max_h = max(ess[:,0])+1
+            for h in range(max_h):
+                hidxs = np.argwhere(ess[:,0] == h)[:,0]
+                output[1].append(ess[hidxs][:,1])
+
+            reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind)
+            reg = pr[reg_ind]
+            max_h = max(reg[:,0])+1
+            for h in range(max_h):
+                hidxs = np.argwhere(reg[:,0] == h)[:,0]
+                output[0].append(reg[hidxs][:,1:])
         else:
             print("cofaces_of_persistence_pairs function requires persistence function"
                   " to be launched first.")
-        return np.array(persistence_result)
+        return output
 
     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 59024212..32e6ee9c 100644
--- a/src/python/include/Persistent_cohomology_interface.h
+++ b/src/python/include/Persistent_cohomology_interface.h
@@ -16,6 +16,7 @@
 #include <vector>
 #include <utility>  // for std::pair
 #include <algorithm>  // for sort
+#include <unordered_map>
 
 namespace Gudhi {
 
@@ -80,8 +81,9 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol
     }
     // 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));
+    std::unordered_map<int, int> order;  
+    //std::sort(max_splx.begin(), max_splx.end()); 
+    for (unsigned int i = 0; i < max_splx.size(); i++)  order.emplace(max_splx[i], i);
 
     std::vector<std::vector<int>> persistence_pairs;
     for (auto pair : pairs) {
diff --git a/src/python/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py
index dd7653c2..5c59db8f 100755
--- a/src/python/test/test_cubical_complex.py
+++ b/src/python/test/test_cubical_complex.py
@@ -151,4 +151,9 @@ def test_connected_sublevel_sets():
 def test_cubical_generators():
     cub = CubicalComplex(top_dimensional_cells = [[0, 0, 0], [0, 1, 0], [0, 0, 0]])
     cub.persistence()
-    assert np.array_equal(cub.cofaces_of_persistence_pairs(), np.array([[1, 7, 4], [0, 8, -1]]))
+    g = cub.cofaces_of_persistence_pairs()
+    assert len(g[0]) == 2
+    assert len(g[1]) == 1
+    assert np.array_equal(g[0][0], np.empty(shape=[0,2]))
+    assert np.array_equal(g[0][1], np.array([[7, 4]]))
+    assert np.array_equal(g[1][0], np.array([8]))