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;
+  /// If true, assumes that Filtration_value is vector-like instead of float-like. This also assumes that Filtration_values is a class, which has a push_to method to push a filtration value $x$ onto $this>=0$. 
+  static const bool is_multi_parameter;
 };
 

src/Simplex_tree/include/gudhi/Simplex_tree.h

@@ -154,16 +154,18 @@ 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; }
+    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;
@@ -586,7 +588,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();
   }
@@ -614,18 +616,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);
@@ -848,7 +857,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();
@@ -914,7 +923,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);
 
@@ -943,7 +952,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);
 
@@ -972,7 +981,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
@@ -989,8 +998,20 @@ 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)) { 
+        if constexpr (SimplexTreeOptions::is_multi_parameter){
+          // By default, does nothing and assumes that the user is smart.
+          if (filt < filtration(simplex_one)){
+            // placeholder for comparable filtrations
+          }
+          else{
+            // placeholder for incomparable filtrations
+          }
+        }
+        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();
@@ -1335,7 +1356,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));
@@ -1689,7 +1710,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) {
@@ -1901,12 +1922,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),
+      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 = 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)) {
@@ -1941,7 +1972,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);
@@ -1951,7 +1982,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();
     bool modified = false;
     bool emptied = false;
@@ -2383,7 +2414,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.
   }
@@ -2394,7 +2425,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);
@@ -2560,6 +2591,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.
@@ -2611,6 +2654,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
@@ -2628,6 +2672,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
@@ -2645,6 +2690,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
+