2 files changed, 127 insertions, 6 deletions

diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 708cdef9..3ba838a7 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -30,6 +30,7 @@
 
 #include <gudhi/reader_utils.h>
 #include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/Debug_utils.h>
 
 #include <boost/container/flat_map.hpp>
 #include <boost/iterator/transform_iterator.hpp>
@@ -43,10 +44,11 @@
 #include <utility>
 #include <vector>
 #include <functional>  // for greater<>
+#include <stdexcept>
+#include <limits>  // Inf
 #include <initializer_list>
 
 namespace Gudhi {
-
 /** \defgroup simplex_tree Filtered Complexes
  *
  * A simplicial complex \f$\mathbf{K}\f$
@@ -446,6 +448,15 @@ class Simplex_tree {
     }
   }
 
+  /** \brief Sets the filtration value of a simplex.
+   *
+   * No action if called on the null_simplex*/
+  void assign_filtration(Simplex_handle sh, Filtration_value fv) {
+    if (sh != null_simplex()) {
+      sh->second.assign_filtration(fv);
+    }
+  }
+
   /** \brief Returns an upper bound of the filtration values of the simplices. */
   Filtration_value filtration() const {
     return threshold_;
@@ -515,7 +526,7 @@ class Simplex_tree {
     return dimension_;
   }
 
-  /** \brief Returns true iff the node in the simplex tree pointed by
+  /** \brief Returns true if the node in the simplex tree pointed by
    * sh has children.*/
   bool has_children(Simplex_handle sh) const {
     return (sh->second.children()->parent() == sh->first);
@@ -718,7 +729,7 @@ class Simplex_tree {
     } else if (the_simplex.size() == 1) {
       // When reaching the end of recursivity, vector of simplices shall be empty and filled on back recursive
       if ((to_be_inserted.size() != 0) || (to_be_propagated.size() != 0)) {
-        std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Error vector not empty" << std::endl;
+        std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Error vector not empty\n";
         exit(-1);
       }
       std::vector<Vertex_handle> first_simplex(1, the_simplex.back());
@@ -727,7 +738,7 @@ class Simplex_tree {
 
       insert_result = insert_vertex_vector(first_simplex, filtration);
     } else {
-        std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Recursivity error" << std::endl;
+        std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Recursivity error\n";
         exit(-1);
     }
     return insert_result;
@@ -740,7 +751,6 @@ class Simplex_tree {
     sh->second.assign_key(key);
   }
 
- public:
   /** Returns the two Simplex_handle corresponding to the endpoints of
    * and edge. sh must point to a 1-dimensional simplex. This is an
    * optimized version of the boundary computation. */
@@ -1105,6 +1115,114 @@ class Simplex_tree {
       os << filtration(sh) << " \n";
     }
   }
+  
+ public:
+  /** \brief Browse the simplex tree to ensure the filtration is not decreasing.
+   * The simplex tree is browsed starting from the root until the leaf, and the filtration values are set with their
+   * parent value (increased), in case the values are decreasing.
+   * @return The filtration modification information.
+   * \warning Some simplex tree functions require the filtration to be valid. `make_filtration_non_decreasing()`
+   * function is not launching `initialize_filtration()` but returns the filtration modification information. If the
+   * complex has changed , please call `initialize_filtration()` to recompute it.
+   */
+  bool make_filtration_non_decreasing() {
+    bool modified = false;
+    for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) {
+      if (has_children(sh)) {
+        modified |= rec_make_filtration_non_decreasing(sh->second.children(), sh->second.filtration());
+      }
+    }
+    return modified;
+  }
+  
+ private:
+  /** \brief Recursively Browse the simplex tree to ensure the filtration is not decreasing.
+   * @param[in] sib Siblings to be parsed.
+   * @param[in] upper_filtration Upper level filtration value in the simplex tree.
+   * @return The filtration modification information in order to trigger initialize_filtration.
+   */
+  bool rec_make_filtration_non_decreasing(Siblings * sib, Filtration_value upper_filtration) {
+    bool modified = false;
+    for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
+      if (sh->second.filtration() < upper_filtration) {
+        // Store the filtration modification information
+        modified = true;
+        sh->second.assign_filtration(upper_filtration);
+      }
+      if (has_children(sh)) {
+        modified |= rec_make_filtration_non_decreasing(sh->second.children(), sh->second.filtration());
+      }
+    }
+    // Make the modified information to be traced by upper call
+    return modified;
+  }
+
+ public:
+  /** \brief Prune above filtration value given as parameter.
+   * @param[in] filtration Maximum threshold value.
+   * \warning threshold_ is set from filtration given as parameter.
+   * \warning The filtration must be valid. If the filtration has not been initialized yet, the method initializes it
+   * (i.e. order the simplices). If the complex has changed since the last time the filtration was initialized, please
+   * call `initialize_filtration()` to recompute it.
+   */
+  void prune_above_filtration(Filtration_value filtration) {
+std::cout << "prune_above_filtration - filtration=" << filtration << std::endl;
+    // No action if filtration is not stored
+    if (Options::store_filtration) {
+      if (filtration < threshold_) {
+        threshold_ = filtration;
+        // Initialize filtration_vect_ if required
+        if (filtration_vect_.empty()) {
+          initialize_filtration();
+        }
+        
+        std::vector<std::vector<Vertex_handle>> simplex_list_to_removed;
+        // Loop in reverse mode until threshold is reached
+        // Do not erase while looping, because removing is shifting data in a flat_map
+        for (auto f_simplex = filtration_vect_.rbegin();
+             (f_simplex != filtration_vect_.rend()) && ((*f_simplex)->second.filtration() > threshold_);
+             f_simplex++) {
+          std::vector<Vertex_handle> simplex_to_remove;
+          for (auto vertex : simplex_vertex_range(*f_simplex))
+            simplex_to_remove.insert(simplex_to_remove.begin(), vertex);
+          simplex_list_to_removed.push_back(simplex_to_remove);
+        }
+        for (auto simplex_to_remove : simplex_list_to_removed) {
+          Simplex_handle sh = find_simplex(simplex_to_remove);
+          if (sh != null_simplex())
+            remove_maximal_simplex(sh);
+        }
+        // Re-initialize filtration_vect_ if dta were removed, because removing is shifting data in a flat_map
+        if (simplex_list_to_removed.size() > 0)
+          initialize_filtration();
+      }
+    }
+  }
+
+  /** \brief Remove a maximal simplex.
+   * @param[in] sh Simplex handle on the maximal simplex to remove.
+   * \pre Please check the simplex has no coface before removing it.
+   * \warning In debug mode, the exception std::invalid_argument is thrown if sh has children.
+   * \warning Be aware that removing is shifting data in a flat_map (initialize_filtration to be done).
+   */
+  void remove_maximal_simplex(Simplex_handle sh) {
+    // Guarantee the simplex has no children
+    GUDHI_CHECK(has_children(sh),
+                std::invalid_argument ("Simplex_tree::remove_maximal_simplex - argument has children"));
+    
+    // Simplex is a leaf, it means the child is the Siblings owning the leaf
+    Siblings* child = sh->second.children();
+
+    if ((child->size() > 1) || (child == root())) {
+      // Not alone, just remove it from members
+      // Special case when child is the root of the simplex tree, just remove it from members
+      child->erase(sh->first);
+    } else {
+      // Sibling is emptied : must be deleted, and its parent must point on his own Sibling
+      child->oncles()->members().at(child->parent()).assign_children(child->oncles());
+      delete child;
+    }
+  }
 
  private:
   Vertex_handle null_vertex_;
@@ -1120,7 +1238,6 @@ class Simplex_tree {
 };
 
 // Print a Simplex_tree in os.
-
 template<typename...T>
 std::ostream& operator<<(std::ostream & os, Simplex_tree<T...> & st) {
   for (auto sh : st.filtration_simplex_range()) {
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h
index 072afc8d..8942d92c 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h
@@ -116,6 +116,10 @@ class Simplex_tree_siblings {
     return members_.size();
   }
 
+  void erase(const Vertex_handle vh) {
+    members_.erase(vh);
+  }
+
   Simplex_tree_siblings * oncles_;
   Vertex_handle parent_;
   Dictionary members_;