1 files changed, 44 insertions, 0 deletions
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();