Merge pull request #204 from MathieuCarriere/generators

Generators

5 files changed, 185 insertions, 0 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 e6a78a6d..f8f80ded 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
@@ -13,6 +13,8 @@
 
 #include <gudhi/Bitmap_cubical_complex/counter.h>
 
+#include <boost/config.hpp>
+
 #include <iostream>
 #include <vector>
 #include <string>
@@ -110,6 +112,16 @@ class Bitmap_cubical_complex_base {
   virtual inline std::vector<std::size_t> get_coboundary_of_a_cell(std::size_t cell) const;
 
   /**
+   * 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 can be a cell of any dimension (vertex, edge, etc).
+   * On the other hand, the output is always indicating the position of
+   * a top-dimensional cube in the data structure.
+   * \pre The filtration values are assigned as per `impose_lower_star_filtration()`.
+   **/
+  inline size_t get_top_dimensional_coface_of_a_cell(size_t splx);
+
+  /**
   * This procedure compute incidence numbers between cubes. For a cube \f$A\f$ of
   * dimension n and a cube \f$B \subset A\f$ of dimension n-1, an incidence
   * between \f$A\f$ and \f$B\f$ is the integer with which \f$B\f$ appears in the boundary of \f$A\f$.
@@ -603,6 +615,19 @@ void Bitmap_cubical_complex_base<T>::setup_bitmap_based_on_top_dimensional_cells
 }
 
 template <typename T>
+size_t Bitmap_cubical_complex_base<T>::get_top_dimensional_coface_of_a_cell(size_t splx) {
+  if (this->get_dimension_of_a_cell(splx) == this->dimension()){return splx;}
+  else{
+    for (auto v : this->get_coboundary_of_a_cell(splx)){
+      if(this->get_cell_data(v) == this->get_cell_data(splx)){
+        return this->get_top_dimensional_coface_of_a_cell(v);
+      }
+    }
+  }
+  BOOST_UNREACHABLE_RETURN(-2);
+}
+
+template <typename T>
 Bitmap_cubical_complex_base<T>::Bitmap_cubical_complex_base(const std::vector<unsigned>& sizes_in_following_directions,
                                                             const std::vector<T>& top_dimensional_cells) {
   this->setup_bitmap_based_on_top_dimensional_cells_list(sizes_in_following_directions, top_dimensional_cells);
diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx
index 007abcb6..ca979eda 100644
--- a/src/python/gudhi/cubical_complex.pyx
+++ b/src/python/gudhi/cubical_complex.pyx
@@ -37,6 +37,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
         Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max)
         void compute_persistence(int homology_coeff_field, double min_persistence)
         vector[pair[int, pair[double, double]]] get_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)
@@ -170,6 +171,58 @@ cdef class CubicalComplex:
         self.compute_persistence(homology_coeff_field, min_persistence)
         return self.pcohptr.get_persistence()
 
+    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,
+            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).
+        """
+
+        assert self.pcohptr != NULL, "compute_persistence() must be called before cofaces_of_persistence_pairs()"
+
+        cdef vector[vector[int]] persistence_result
+        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:])
+
+        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 246a3a02..06309772 100644
--- a/src/python/gudhi/periodic_cubical_complex.pyx
+++ b/src/python/gudhi/periodic_cubical_complex.pyx
@@ -34,6 +34,7 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
         Periodic_cubical_complex_persistence_interface(Periodic_cubical_complex_base_interface * st, bool persistence_dim_max)
         void compute_persistence(int homology_coeff_field, double min_persistence)
         vector[pair[int, pair[double, double]]] get_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)
@@ -175,6 +176,58 @@ cdef class PeriodicCubicalComplex:
         self.compute_persistence(homology_coeff_field, min_persistence)
         return self.pcohptr.get_persistence()
 
+    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,
+            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 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 0de9bd5c..e5a3dfba 100644
--- a/src/python/include/Persistent_cohomology_interface.h
+++ b/src/python/include/Persistent_cohomology_interface.h
@@ -13,9 +13,11 @@
 
 #include <gudhi/Persistent_cohomology.h>
 
+#include <cstdlib>
 #include <vector>
 #include <utility>  // for std::pair
 #include <algorithm>  // for sort
+#include <unordered_map>
 
 namespace Gudhi {
 
@@ -67,6 +69,48 @@ persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomol
     return persistence;
   }
 
+  // This function computes the top-dimensional cofaces associated to the positive and negative
+  // simplices of a cubical complex. The output format is a vector of vectors of three integers,
+  // which are [homological dimension, index of top-dimensional coface of positive simplex,
+  // index of top-dimensional coface of negative simplex]. If the topological feature is essential,
+  // then the index of top-dimensional coface of negative simplex is arbitrarily set to -1.
+  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<std::size_t> 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::unordered_map<std::size_t, 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) {
+      int h = stptr_->dimension(get<0>(pair));
+      // Recursively get the top-dimensional cell / coface associated to the persistence generator
+      std::size_t face0 = stptr_->get_top_dimensional_coface_of_a_cell(get<0>(pair));
+      // Retrieve the index of the corresponding top-dimensional cell in the input data
+      int splx0 = order[face0];
+
+      int splx1 = -1;
+      if (get<1>(pair) != stptr_->null_simplex()){
+        // Recursively get the top-dimensional cell / coface associated to the persistence generator
+        std::size_t face1 = stptr_->get_top_dimensional_coface_of_a_cell(get<1>(pair));
+        // Retrieve the index of the corresponding top-dimensional cell in the input data
+        splx1 = order[face1];
+      }
+      persistence_pairs.push_back({ h, splx0, splx1 });
+    }
+    return persistence_pairs;
+  }
+
   std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs() {
     std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs;
     auto const& pairs = Base::get_persistent_pairs();
diff --git a/src/python/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py
index fce4875c..5c59db8f 100755
--- a/src/python/test/test_cubical_complex.py
+++ b/src/python/test/test_cubical_complex.py
@@ -147,3 +147,13 @@ def test_connected_sublevel_sets():
                                  periodic_dimensions = periodic_dimensions)
     assert cub.persistence() == [(0, (2.0, float("inf")))]
     assert cub.betti_numbers() == [1, 0, 0]
+
+def test_cubical_generators():
+    cub = CubicalComplex(top_dimensional_cells = [[0, 0, 0], [0, 1, 0], [0, 0, 0]])
+    cub.persistence()
+    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]))