Multiparameter simplextree

src/Persistent_cohomology/test/betti_numbers_unit_test.cpp

@@ -19,6 +19,8 @@ struct MiniSTOptions : Gudhi::Simplex_tree_options_full_featured {
   static const bool store_key = true;
   // I have few vertices
   typedef short Vertex_handle;
+
+  static const bool is_multi_parameter = false;
 };
 
 using Mini_simplex_tree = Gudhi::Simplex_tree<MiniSTOptions>;

src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp

@@ -177,6 +177,7 @@ struct MiniSTOptions {
   static const bool contiguous_vertices = false;
   static const bool link_nodes_by_label = false;
   static const bool stable_simplex_handles = false;
+  static const bool is_multi_parameter = false;
 };
 
 using Mini_simplex_tree = Gudhi::Simplex_tree<MiniSTOptions>;

src/Simplex_tree/benchmark/simplex_tree_cofaces_benchmark.cpp

@@ -91,6 +91,7 @@ struct Simplex_tree_options_stable_simplex_handles {
   static const bool contiguous_vertices = false;
   static const bool link_nodes_by_label = true;
   static const bool stable_simplex_handles = true;
+  static const bool is_multi_parameter = false;
 };
 
 int main(int argc, char *argv[]) {

src/Simplex_tree/concept/SimplexTreeOptions.h

@@ -31,5 +31,7 @@ struct SimplexTreeOptions {
   static const bool link_nodes_by_label;
   /// If true, Simplex_handle will not be invalidated after insertions or removals.
   static const bool stable_simplex_handles;
+
+  static const bool is_multi_parameter;
 };
 

src/Simplex_tree/include/gudhi/Simplex_tree.h

@@ -151,16 +151,19 @@ class Simplex_tree {
       Key_simplex_base;
 
   struct Filtration_simplex_base_real {
-    Filtration_simplex_base_real() : filt_(0) {}
-    void assign_filtration(Filtration_value f) { filt_ = f; }
-    Filtration_value filtration() const { return filt_; }
+    Filtration_simplex_base_real() : filt_{} {}
+    void assign_filtration(const Filtration_value& f) { filt_ = f; }
+    // Filtration_value filtration() const { return filt_; }
+    const Filtration_value& filtration() const { return filt_; }
+    Filtration_value& filtration() { return filt_; }
    private:
     Filtration_value filt_;
   };
   struct Filtration_simplex_base_dummy {
     Filtration_simplex_base_dummy() {}
     void assign_filtration(Filtration_value GUDHI_CHECK_code(f)) { GUDHI_CHECK(f == 0, "filtration value specified for a complex that does not store them"); }
-    Filtration_value filtration() const { return 0; }
+    const Filtration_value& filtration() const  { return null_value; }
+    static constexpr Filtration_value null_value={};
   };
   typedef typename std::conditional<Options::store_filtration, Filtration_simplex_base_real,
     Filtration_simplex_base_dummy>::type Filtration_simplex_base;
@@ -576,7 +579,7 @@ class Simplex_tree {
    *
    * Same as `filtration()`, but does not handle `null_simplex()`.
    */
-  static Filtration_value filtration_(Simplex_handle sh) {
+  static const Filtration_value& filtration_(Simplex_handle sh) {
     GUDHI_CHECK (sh != null_simplex(), "null simplex");
     return sh->second.filtration();
   }
@@ -604,18 +607,25 @@ class Simplex_tree {
    * Called on the null_simplex, it returns infinity.
    * If SimplexTreeOptions::store_filtration is false, returns 0.
    */
-  static Filtration_value filtration(Simplex_handle sh) {
+  static const Filtration_value& filtration(Simplex_handle sh){
     if (sh != null_simplex()) {
       return sh->second.filtration();
     } else {
-      return std::numeric_limits<Filtration_value>::infinity();
+      return inf_;
+    }
+  }
+  static Filtration_value& filtration_mutable(Simplex_handle sh){
+    if (sh != null_simplex()) {
+      return sh->second.filtration();
+    } else {
+      return inf_;
     }
   }
 
   /** \brief Sets the filtration value of a simplex.
    * \exception std::invalid_argument In debug mode, if sh is a null_simplex.
    */
-  void assign_filtration(Simplex_handle sh, Filtration_value fv) {
+  void assign_filtration(Simplex_handle sh, const Filtration_value& fv) {
     GUDHI_CHECK(sh != null_simplex(),
                 std::invalid_argument("Simplex_tree::assign_filtration - cannot assign filtration on null_simplex"));
     sh->second.assign_filtration(fv);
@@ -829,7 +839,7 @@ class Simplex_tree {
   */
   template <class RandomVertexHandleRange = std::initializer_list<Vertex_handle>>
   std::pair<Simplex_handle, bool> insert_simplex_raw(const RandomVertexHandleRange& simplex,
-                                                     Filtration_value filtration) {
+                                                     const Filtration_value& filtration) {
     Siblings * curr_sib = &root_;
     std::pair<Simplex_handle, bool> res_insert;
     auto vi = simplex.begin();
@@ -895,7 +905,7 @@ class Simplex_tree {
    * .end() return input iterators, with 'value_type' Vertex_handle. */
   template<class InputVertexRange = std::initializer_list<Vertex_handle>>
   std::pair<Simplex_handle, bool> insert_simplex(const InputVertexRange & simplex,
-                                                 Filtration_value filtration = 0) {
+                                                 const Filtration_value& filtration = {}) {
     auto first = std::begin(simplex);
     auto last = std::end(simplex);
 
@@ -924,7 +934,7 @@ class Simplex_tree {
    */
   template<class InputVertexRange = std::initializer_list<Vertex_handle>>
   std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(const InputVertexRange& Nsimplex,
-                    Filtration_value filtration = 0) {
+				    const Filtration_value& filtration = {}) {
     auto first = std::begin(Nsimplex);
     auto last = std::end(Nsimplex);
 
@@ -953,7 +963,7 @@ class Simplex_tree {
   std::pair<Simplex_handle, bool> rec_insert_simplex_and_subfaces_sorted(Siblings* sib,
                                                                            ForwardVertexIterator first,
                                                                            ForwardVertexIterator last,
-                                                                           Filtration_value filt) {
+  	                                                                     const Filtration_value& filt) {
     // An alternative strategy would be:
     // - try to find the complete simplex, if found (and low filtration) exit
     // - insert all the vertices at once in sib
@@ -970,8 +980,28 @@ class Simplex_tree {
     Simplex_handle simplex_one = insertion_result.first;
     bool one_is_new = insertion_result.second;
     if (!one_is_new) {
-      if (filtration(simplex_one) > filt) {
-        assign_filtration(simplex_one, filt);
+      if (!(filtration(simplex_one) <= filt)) { // TODO : For multipersistence, it's not clear what should be the default, especially for multicritical filtrations
+        if constexpr (SimplexTreeOptions::is_multi_parameter){
+          if (filt < filtration(simplex_one)){
+          //   assign_filtration(simplex_one, filt);
+          // I don't really like this behavior. 
+            // It prevents inserting simplices by default at the smallest possible position after its childrens.
+            // e.g., if (python) we type : st.insert([0], [1,0]), st.insert([1], [0,1]), st.insert([0,1])
+            // we may want st.filtration([0,1]) to be [1,1]. (maybe after a make_filtration_decreasing from user)
+            // Furthermore, this may more sense as default value -> 0 can erase filtration values of childrens... 
+          }
+          else{ // As multicritical filtrations are not well supported yet, its not safe to concatenate yet.
+            // two incomparables filtrations values -> we concatenate 
+            // std::cout << "incomparable -> concatenate" << std::endl;
+            // Filtration_value new_filtration = filtration(simplex_one);
+            // new_filtration.insert_new(filt); // we assume then that Filtration_value has a insert_new method.
+            // assign_filtration(simplex_one, new_filtration);
+          }
+        }
+        else{ // non-multiparameter
+          assign_filtration(simplex_one, filt);
+        }
+
       } else {
         // FIXME: this interface makes no sense, and it doesn't seem to be tested.
         insertion_result.first = null_simplex();
@@ -1316,7 +1346,7 @@ class Simplex_tree {
    * The complex does not need to be empty before calling this function. However, if a vertex is
    * already present, its filtration value is not modified, unlike with other insertion functions. */
   template <class VertexRange>
-  void insert_batch_vertices(VertexRange const& vertices, Filtration_value filt = 0) {
+  void insert_batch_vertices(VertexRange const& vertices, const Filtration_value& filt ={}) {
     auto verts = vertices | boost::adaptors::transformed([&](auto v){
         return Dit_value_t(v, Node(&root_, filt)); });
     root_.members_.insert(boost::begin(verts), boost::end(verts));
@@ -1396,7 +1426,7 @@ class Simplex_tree {
   static void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection,
                            Dictionary_it begin1, Dictionary_it end1,
                            Dictionary_it begin2, Dictionary_it end2,
-                           Filtration_value filtration_) {
+                           const Filtration_value& filtration_) {
     if (begin1 == end1 || begin2 == end2)
       return;  // ----->>
     while (true) {
@@ -1603,12 +1633,22 @@ class Simplex_tree {
       if (dim == 0) return;
       // Find the maximum filtration value in the border
       Boundary_simplex_range&& boundary = boundary_simplex_range(sh);
-      Boundary_simplex_iterator max_border = std::max_element(std::begin(boundary), std::end(boundary),
+      typename SimplexTreeOptions::Filtration_value max_filt_border_value;
+      if constexpr (SimplexTreeOptions::is_multi_parameter){
+        // in that case, we assume that Filtration_value has a `push_to` member to handle this.
+        max_filt_border_value = typename SimplexTreeOptions::Filtration_value(this->number_of_parameters_);
+        for (auto &face_sh : boundary){
+          max_filt_border_value.push_to(filtration(face_sh)); // pushes the value of max_filt_border_value to reach simplex' filtration
+        }
+      }
+      else{
+        Boundary_simplex_iterator max_border = std::max_element(std::begin(boundary), std::end(boundary),
                                                               [](Simplex_handle sh1, Simplex_handle sh2) {
                                                                 return filtration(sh1) < filtration(sh2);
                                                               });
-
-      Filtration_value max_filt_border_value = filtration(*max_border);
+        max_filt_border_value = filtration(*max_border);
+      }
+
       // Replacing if(f<max) with if(!(f>=max)) would mean that if f is NaN, we replace it with the max of the children.
       // That seems more useful than keeping NaN.
       if (!(sh->second.filtration() >= max_filt_border_value)) {
@@ -1643,7 +1683,7 @@ class Simplex_tree {
    * than it was before. However, `upper_bound_dimension()` will return the old value, which remains a valid upper
    * bound. If you care, you can call `dimension()` to recompute the exact dimension.
    */
-  bool prune_above_filtration(Filtration_value filtration) {
+  bool prune_above_filtration(const Filtration_value& filtration) {
     if (std::numeric_limits<Filtration_value>::has_infinity && filtration == std::numeric_limits<Filtration_value>::infinity())
       return false;  // ---->>
     bool modified = rec_prune_above_filtration(root(), filtration);
@@ -1653,7 +1693,7 @@ class Simplex_tree {
   }
 
  private:
-  bool rec_prune_above_filtration(Siblings* sib, Filtration_value filt) {
+  bool rec_prune_above_filtration(Siblings* sib, const Filtration_value& filt) {
     auto&& list = sib->members();
     auto last = std::remove_if(list.begin(), list.end(), [this,filt](Dit_value_t& simplex) {
         if (simplex.second.filtration() <= filt) return false;
@@ -2063,7 +2103,7 @@ class Simplex_tree {
    * @param[in] filt_value The new filtration value.
    * @param[in] min_dim The minimal dimension. Default value is 0.
    */
-  void reset_filtration(Filtration_value filt_value, int min_dim = 0) {
+  void reset_filtration(const Filtration_value& filt_value, int min_dim = 0) {
     rec_reset_filtration(&root_, filt_value, min_dim);
     clear_filtration(); // Drop the cache.
   }
@@ -2074,7 +2114,7 @@ class Simplex_tree {
    * @param[in] filt_value The new filtration value.
    * @param[in] min_depth The minimal depth.
    */
-  void rec_reset_filtration(Siblings * sib, Filtration_value filt_value, int min_depth) {
+  void rec_reset_filtration(Siblings * sib, const Filtration_value& filt_value, int min_depth) {
     for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
       if (min_depth <= 0) {
         sh->second.assign_filtration(filt_value);
@@ -2239,6 +2279,18 @@ class Simplex_tree {
   /** \brief Upper bound on the dimension of the simplicial complex.*/
   int dimension_;
   bool dimension_to_be_lowered_ = false;
+
+//MULTIPERS STUFF
+public: 
+  void set_number_of_parameters(int num){
+		number_of_parameters_ = num;
+	}
+	int get_number_of_parameters() const{
+		return number_of_parameters_;
+	}
+  inline static Filtration_value inf_ = std::numeric_limits<Filtration_value>::infinity();
+private:
+	int number_of_parameters_;
 };
 
 // Print a Simplex_tree in os.
@@ -2290,6 +2342,7 @@ struct Simplex_tree_options_full_featured {
   static const bool contiguous_vertices = false;
   static const bool link_nodes_by_label = false;
   static const bool stable_simplex_handles = false;
+  static const bool is_multi_parameter = false;
 };
 
 /** Model of SimplexTreeOptions, faster than `Simplex_tree_options_full_featured` but note the unsafe
@@ -2307,6 +2360,7 @@ struct Simplex_tree_options_fast_persistence {
   static const bool contiguous_vertices = true;
   static const bool link_nodes_by_label = false;
   static const bool stable_simplex_handles = false;
+  static const bool is_multi_parameter = false;
 };
 
 /** Model of SimplexTreeOptions, faster cofaces than `Simplex_tree_options_full_featured`, note the
@@ -2324,6 +2378,8 @@ struct Simplex_tree_options_fast_cofaces {
   static const bool contiguous_vertices = false;
   static const bool link_nodes_by_label = true;
   static const bool stable_simplex_handles = false;
+  static const bool is_multi_parameter = false;
+
 };
 
 /** @}*/  // end addtogroup simplex_tree

src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h

@@ -41,7 +41,7 @@ struct GUDHI_EMPTY_BASE_CLASS_OPTIMIZATION Simplex_tree_node_explicit_storage :
   typedef typename SimplexTree::Simplex_key Simplex_key;
 
   Simplex_tree_node_explicit_storage(Siblings * sib = nullptr,
-                                     Filtration_value filtration = 0)
+                                     const Filtration_value& filtration = {}) 
       : children_(sib) {
     this->assign_filtration(filtration);
   }

src/Simplex_tree/test/simplex_tree_ctor_and_move_unit_test.cpp

@@ -36,6 +36,7 @@ struct Simplex_tree_options_stable_simplex_handles {
   static const bool contiguous_vertices = false;
   static const bool link_nodes_by_label = true;
   static const bool stable_simplex_handles = true;
+  static const bool is_multi_parameter = false;
 };
 
 typedef boost::mpl::list<Simplex_tree<>,

src/Simplex_tree/test/simplex_tree_graph_expansion_unit_test.cpp

@@ -31,6 +31,7 @@ struct Simplex_tree_options_stable_simplex_handles {
   static const bool contiguous_vertices = false;
   static const bool link_nodes_by_label = false;
   static const bool stable_simplex_handles = true;
+  static const bool is_multi_parameter = false;
 };
 
 typedef boost::mpl::list<Simplex_tree<>,

src/Simplex_tree/test/simplex_tree_make_filtration_non_decreasing_unit_test.cpp

@@ -33,6 +33,7 @@ struct Simplex_tree_options_stable_simplex_handles {
   static const bool contiguous_vertices = false;
   static const bool link_nodes_by_label = false;
   static const bool stable_simplex_handles = true;
+  static const bool is_multi_parameter = false;
 };
 
 typedef boost::mpl::list<Simplex_tree<>,

src/Simplex_tree/test/simplex_tree_unit_test.cpp

@@ -41,6 +41,7 @@ struct Simplex_tree_options_stable_simplex_handles {
   static const bool contiguous_vertices = false;
   static const bool link_nodes_by_label = true;
   static const bool stable_simplex_handles = true;
+  static const bool is_multi_parameter = false;
 };
 
 typedef boost::mpl::list<Simplex_tree<>,

src/python/CMakeLists.txt

@@ -55,6 +55,7 @@ if(PYTHONINTERP_FOUND)
     # Cython modules
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'off_utils', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'simplex_tree', ")
+    set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'simplex_tree_multi', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'edge_collapse', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'rips_complex', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'cubical_complex', ")
@@ -165,6 +166,7 @@ if(PYTHONINTERP_FOUND)
 
     set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_utils', ")
     set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'simplex_tree', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'simplex_tree_multi', ")
     set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'rips_complex', ")
     set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'cubical_complex', ")
     set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'periodic_cubical_complex', ")

src/python/gudhi/multiparameter_edge_collapse.py

@@ -0,0 +1,29 @@
+from tqdm import tqdm
+from filtration_domination import remove_strongly_filtration_dominated, remove_filtration_dominated
+
+def _collapse_edge_list(edges, num:int=0, full:bool=False, strong:bool=False, progress:bool=False):
+	"""
+	Given an edge list defining a 1 critical 2 parameter 1 dimensional simplicial complex, simplificates this filtered simplicial complex, using filtration-domination's edge collapser.
+	"""
+	n = len(edges)
+	if full:
+		num = 100
+	with tqdm(range(num), total=num, desc="Removing edges", disable=not(progress)) as I:
+		for i in I:
+			if strong:
+				edges = remove_strongly_filtration_dominated(edges) # nogil ?
+			else:
+				edges = remove_filtration_dominated(edges)
+			# Prevents doing useless collapses
+			if len(edges) >= n:
+				if full and strong:
+					strong = False
+					n = len(edges)
+					# n = edges.size() # len(edges)
+				else : 
+					break
+			else:
+				n = len(edges)
+				# n = edges.size()
+	return edges
+