Multiparameter persistence, C++ part.

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/example/CMakeLists.txt

@@ -21,6 +21,13 @@ if(TARGET TBB::tbb)
 endif()
 add_test(NAME Simplex_tree_example_mini_simplex_tree COMMAND $<TARGET_FILE:Simplex_tree_example_mini_simplex_tree>)
 
+add_executable ( Simplex_tree_multi_example simplex_tree_multi.cpp )
+if(TARGET TBB::tbb)
+  target_link_libraries(Simplex_tree_multi_example TBB::tbb)
+endif()
+add_test(NAME Simplex_tree_multi_example COMMAND $<TARGET_FILE:Simplex_tree_multi_example>)
+
+
 # An example with Simplex-tree using CGAL alpha_shapes_3
 if(GMP_FOUND AND NOT CGAL_VERSION VERSION_LESS 4.11.0)
   add_executable ( Simplex_tree_example_alpha_shapes_3_from_off example_alpha_shapes_3_simplex_tree_from_off_file.cpp )

src/Simplex_tree/example/simplex_tree_multi.cpp

@@ -0,0 +1,56 @@
+/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ *    Author(s):       David Loiseaux
+ *
+ *    Copyright (C) 2023 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Simplex_tree/Simplex_tree_multi.h>
+
+#include <iostream>
+#include <initializer_list>
+
+struct ST_MULTI {
+public:
+	typedef Gudhi::linear_indexing_tag Indexing_tag;
+	typedef int Vertex_handle;
+	typedef float value_type;
+	using Filtration_value = Gudhi::multiparameter::multi_filtrations::Finitely_critical_multi_filtration<value_type>;
+	typedef std::uint32_t Simplex_key;
+	static const bool store_key = true;
+	static const bool store_filtration = true;
+	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 = true;
+};
+
+using ST = Gudhi::Simplex_tree<ST_MULTI>;
+
+
+int main() {
+  ST st;
+
+  /* Complex to build. */
+  /*    1              */
+  /*    o              */
+  /*   /X\             */
+  /*  o---o---o        */
+  /*  2   0   3        */
+
+  auto triangle012 = {0, 1, 2};
+  auto edge03 = {0, 3};
+  st.insert_simplex_and_subfaces(triangle012, {1,2,3}); // {1,2,3} can be any array-like vector-like
+  st.insert_simplex_and_subfaces(edge03, {4,5,6});
+
+  auto edge02 = {0, 2};
+  ST::Simplex_handle e = st.find(edge02);
+  // Finitely_critical_multi_filtration has an operator<<
+  std::cout << st.filtration(e) << std::endl;
+  assert(st.filtration(st.find(edge03)) == std::vector<float>({4,5,6}));
+
+}

src/Simplex_tree/include/gudhi/Simplex_tree.h

@@ -151,16 +151,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;
@@ -576,7 +578,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 +606,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);
@@ -822,14 +831,16 @@ class Simplex_tree {
    * to the new simplex.
    * If the insertion fails (the simplex is already there), the bool is set to false. If the insertion
    * fails and the simplex already in the complex has a filtration value strictly bigger than 'filtration',
+   * and the simplex tree is not multi parameter (`SimplexTreeOptions::is_multi_parameter == false`),
    * we assign this simplex with the new value 'filtration', and set the Simplex_handle field of the
-   * output pair to the Simplex_handle of the simplex. Otherwise, we set the Simplex_handle part to
-   * null_simplex.
+   * output pair to the Simplex_handle of the simplex. When the simplex tree is multi parameter,
+   * the existing filtration values are not updated. If the insertion fails for other reasons, 
+   * we set the Simplex_handle part to `null_simplex`.
    * 
   */
   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 +906,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 +935,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 +964,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
@@ -969,14 +980,17 @@ 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);
-      } else {
-        // FIXME: this interface makes no sense, and it doesn't seem to be tested.
-        insertion_result.first = null_simplex();
+    if constexpr (!SimplexTreeOptions::is_multi_parameter){ // Ignores the assign part for multiparameter filtrations.
+      if (!one_is_new) {
+        if (filtration(simplex_one) > filt){ 
+            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();
+        }
       }
     }
+
     if (++first == last) return insertion_result;
     if (!has_children(simplex_one))
       // TODO: have special code here, we know we are building the whole subtree from scratch.
@@ -1316,7 +1330,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));
@@ -1403,7 +1417,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) {
@@ -1610,12 +1624,21 @@ 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),
-                                                              [](Simplex_handle sh1, Simplex_handle sh2) {
-                                                                return filtration(sh1) < filtration(sh2);
-                                                              });
+      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); });
+        max_filt_border_value = filtration(*max_border);
+      }
 
-      Filtration_value 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)) {
@@ -1650,7 +1673,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);
@@ -1660,7 +1683,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;
@@ -2070,7 +2093,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.
   }
@@ -2081,7 +2104,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);
@@ -2246,6 +2269,24 @@ class Simplex_tree {
   /** \brief Upper bound on the dimension of the simplicial complex.*/
   int dimension_;
   bool dimension_to_be_lowered_ = false;
+
+  // MULTIPERS STUFF
+ public:
+  /**
+   * \brief Sets the number of parameters of the filtrations if SimplexTreeOptions::is_multi_parameter. 
+   * */
+  void set_number_of_parameters(int num) { number_of_parameters_ = num; }
+  /**
+   * \brief Gets the number of parameters of the filtrations if SimplexTreeOptions::is_multi_parameter. 
+   * */
+  int get_number_of_parameters() const { return number_of_parameters_; }
+
+  inline static Filtration_value inf_ = std::numeric_limits<Filtration_value>::has_infinity() ? 
+      std::numeric_limits<Filtration_value>::infinity() 
+    : std::numeric_limits<Filtration_value>::max(); /**< Default infinite value. */
+
+ private:
+  int number_of_parameters_; /**< Number of parameters of the multi-filtrations when SimplexTreeOptions::is_multi_parameter.-*/
 };
 
 // Print a Simplex_tree in os.
@@ -2297,6 +2338,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
@@ -2314,6 +2356,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
@@ -2331,6 +2374,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