Merge pull request #242 from mglisse/genrips

Generators for Simplex_tree

3 files changed, 112 insertions, 12 deletions

diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 7b39a500..b7fb9002 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -24,6 +24,7 @@
 #include <boost/iterator/transform_iterator.hpp>
 #include <boost/graph/adjacency_list.hpp>
 #include <boost/range/adaptor/reversed.hpp>
+#include <boost/container/static_vector.hpp>
 
 #ifdef GUDHI_USE_TBB
 #include <tbb/parallel_sort.h>
@@ -246,8 +247,8 @@ class Simplex_tree {
    * which is consequenlty
    * equal to \f$(-1)^{\text{dim} \sigma}\f$ the canonical orientation on the simplex.
    */
-  Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) {
-    assert(sh != null_simplex());  // Empty simplex
+  Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) const {
+    GUDHI_CHECK(sh != null_simplex(), "empty simplex");
     return Simplex_vertex_range(Simplex_vertex_iterator(this, sh),
                                 Simplex_vertex_iterator(this));
   }
@@ -450,6 +451,15 @@ class Simplex_tree {
     return true;
   }
 
+  /** \brief Returns the filtration value of a simplex.
+   *
+   * Same as `filtration()`, but does not handle `null_simplex()`.
+   */
+  static Filtration_value filtration_(Simplex_handle sh) {
+    GUDHI_CHECK (sh != null_simplex(), "null simplex");
+    return sh->second.filtration();
+  }
+
  public:
   /** \brief Returns the key associated to a simplex.
    *
@@ -827,7 +837,7 @@ class Simplex_tree {
 
   /** Returns the Siblings containing a simplex.*/
   template<class SimplexHandle>
-  Siblings* self_siblings(SimplexHandle sh) {
+  static Siblings* self_siblings(SimplexHandle sh) {
     if (sh->second.children()->parent() == sh->first)
       return sh->second.children()->oncles();
     else
@@ -1469,6 +1479,66 @@ class Simplex_tree {
     }
   }
 
+  /** \brief Returns a vertex of `sh` that has the same filtration value as `sh` if it exists, and `null_vertex()` otherwise.
+   *
+   * For a lower-star filtration built with `make_filtration_non_decreasing()`, this is a way to invert the process and find out which vertex had its filtration value propagated to `sh`.
+   * If several vertices have the same filtration value, the one it returns is arbitrary. */
+  Vertex_handle vertex_with_same_filtration(Simplex_handle sh) {
+    auto filt = filtration_(sh);
+    for(auto v : simplex_vertex_range(sh))
+      if(filtration_(find_vertex(v)) == filt)
+        return v;
+    return null_vertex();
+  }
+
+  /** \brief Returns an edge of `sh` that has the same filtration value as `sh` if it exists, and `null_simplex()` otherwise.
+   *
+   * For a flag-complex built with `expansion()`, this is a way to invert the process and find out which edge had its filtration value propagated to `sh`.
+   * If several edges have the same filtration value, the one it returns is arbitrary.
+   *
+   * \pre `sh` must have dimension at least 1. */
+  Simplex_handle edge_with_same_filtration(Simplex_handle sh) {
+    // See issue #251 for potential speed improvements.
+    auto&& vertices = simplex_vertex_range(sh); // vertices in decreasing order
+    auto end = std::end(vertices);
+    auto vi = std::begin(vertices);
+    GUDHI_CHECK(vi != end, "empty simplex");
+    auto v0 = *vi;
+    ++vi;
+    GUDHI_CHECK(vi != end, "simplex of dimension 0");
+    if(std::next(vi) == end) return sh; // shortcut for dimension 1
+    boost::container::static_vector<Vertex_handle, 40> suffix;
+    suffix.push_back(v0);
+    auto filt = filtration_(sh);
+    do
+    {
+      Vertex_handle v = *vi;
+      auto&& children1 = find_vertex(v)->second.children()->members_;
+      for(auto w : suffix){
+        // Can we take advantage of the fact that suffix is ordered?
+        Simplex_handle s = children1.find(w);
+        if(filtration_(s) == filt)
+          return s;
+      }
+      suffix.push_back(v);
+    }
+    while(++vi != end);
+    return null_simplex();
+  }
+
+  /** \brief Returns a minimal face of `sh` that has the same filtration value as `sh`.
+   *
+   * For a filtration built with `make_filtration_non_decreasing()`, this is a way to invert the process and find out which simplex had its filtration value propagated to `sh`.
+   * If several minimal (for inclusion) simplices have the same filtration value, the one it returns is arbitrary, and it is not guaranteed to be the one with smallest dimension. */
+  Simplex_handle minimal_simplex_with_same_filtration(Simplex_handle sh) {
+    auto filt = filtration_(sh);
+    // Naive implementation, it can be sped up.
+    for(auto b : boundary_simplex_range(sh))
+      if(filtration_(b) == filt)
+        return minimal_simplex_with_same_filtration(b);
+    return sh; // None of its faces has the same filtration.
+  }
+
  private:
   Vertex_handle null_vertex_;
   /** \brief Total number of simplices in the complex, without the empty simplex.*/
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h
index efccf2f2..9007b6bd 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h
@@ -15,9 +15,7 @@
 
 #include <boost/iterator/iterator_facade.hpp>
 #include <boost/version.hpp>
-#if BOOST_VERSION >= 105600
-# include <boost/container/static_vector.hpp>
-#endif
+#include <boost/container/static_vector.hpp>
 
 #include <vector>
 
@@ -42,13 +40,13 @@ class Simplex_tree_simplex_vertex_iterator : public boost::iterator_facade<
   typedef typename SimplexTree::Siblings Siblings;
   typedef typename SimplexTree::Vertex_handle Vertex_handle;
 
-  explicit Simplex_tree_simplex_vertex_iterator(SimplexTree * st)
+  explicit Simplex_tree_simplex_vertex_iterator(SimplexTree const* st)
       :  // any end() iterator
         sib_(nullptr),
         v_(st->null_vertex()) {
   }
 
-  Simplex_tree_simplex_vertex_iterator(SimplexTree * st, Simplex_handle sh)
+  Simplex_tree_simplex_vertex_iterator(SimplexTree const* st, Simplex_handle sh)
       : sib_(st->self_siblings(sh)),
         v_(sh->first) {
   }
@@ -166,15 +164,11 @@ class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade<
   // Most of the storage should be moved to the range, iterators should be light.
   Vertex_handle last_;  // last vertex of the simplex
   Vertex_handle next_;  // next vertex to push in suffix_
-#if BOOST_VERSION >= 105600
   // 40 seems a conservative bound on the dimension of a Simplex_tree for now,
   // as it would not fit on the biggest hard-drive.
   boost::container::static_vector<Vertex_handle, 40> suffix_;
   // static_vector still has some overhead compared to a trivial hand-made
   // version using std::aligned_storage, or compared to making suffix_ static.
-#else
-  std::vector<Vertex_handle> suffix_;
-#endif
   Siblings * sib_;  // where the next search will start from
   Simplex_handle sh_;  // current Simplex_handle in the boundary
   SimplexTree * st_;  // simplex containing the simplicial complex
diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
index e739ad0a..3872c4ca 100644
--- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp
+++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
@@ -902,5 +902,41 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(insert_duplicated_vertices, typeST, list_of_tested
             << " - num_simplices = " << st.num_simplices() << std::endl;
   BOOST_CHECK(st.dimension() == 1);
   BOOST_CHECK(st.num_simplices() == st.num_vertices() + 1);
+}
 
+BOOST_AUTO_TEST_CASE_TEMPLATE(generators, typeST, list_of_tested_variants) {
+  std::cout << "********************************************************************" << std::endl;
+  std::cout << "TEST FIND GENERATORS" << std::endl;
+  {
+    typeST st;
+    st.insert_simplex_and_subfaces({0,1,2,3,4,5,6},0);
+    st.assign_filtration(st.find({0,2,4}), 10);
+    st.assign_filtration(st.find({1,5}), 20);
+    st.assign_filtration(st.find({1,2,4}), 30);
+    st.assign_filtration(st.find({3}), 5);
+    st.make_filtration_non_decreasing();
+    BOOST_CHECK(st.filtration(st.find({1,2}))==0);
+    BOOST_CHECK(st.filtration(st.find({0,1,2,3,4}))==30);
+    BOOST_CHECK(st.minimal_simplex_with_same_filtration(st.find({0,1,2,3,4,5}))==st.find({1,2,4}));
+    BOOST_CHECK(st.minimal_simplex_with_same_filtration(st.find({0,2,3}))==st.find({3}));
+    auto s=st.minimal_simplex_with_same_filtration(st.find({0,2,6}));
+    BOOST_CHECK(s==st.find({0})||s==st.find({2})||s==st.find({6}));
+    BOOST_CHECK(st.vertex_with_same_filtration(st.find({2}))==2);
+    BOOST_CHECK(st.vertex_with_same_filtration(st.find({1,5}))==st.null_vertex());
+    BOOST_CHECK(st.vertex_with_same_filtration(st.find({5,6}))>=5);
+  }
+  {
+    typeST st;
+    st.insert_simplex_and_subfaces({0,1}, 8);
+    st.insert_simplex_and_subfaces({0,2}, 10);
+    st.insert_simplex_and_subfaces({3,4}, 6);
+    st.insert_simplex_and_subfaces({1,2}, 5);
+    st.insert_simplex_and_subfaces({1,5}, 4);
+    st.insert_simplex_and_subfaces({0,5}, 3);
+    st.insert_simplex_and_subfaces({2,5}, 2);
+    st.insert_simplex_and_subfaces({1,3}, 9);
+    st.expansion(50);
+    BOOST_CHECK(st.edge_with_same_filtration(st.find({0,1,2,5}))==st.find({0,2}));
+    BOOST_CHECK(st.edge_with_same_filtration(st.find({1,5}))==st.find({1,5}));
+  }
 }