update output

1 files changed, 28 insertions, 5 deletions

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.