1 files changed, 60 insertions, 44 deletions

diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
index 17a237a7..94a125c1 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
@@ -45,9 +45,7 @@ bool Skeleton_blocker_complex<SkeletonBlockerDS>::is_popable_blocker(Blocker_han
   assert(this->contains_blocker(*sigma));
   Skeleton_blocker_link_complex<Skeleton_blocker_complex> link_blocker_sigma;
   build_link_of_blocker(*this, *sigma, link_blocker_sigma);
-
-  bool res = link_blocker_sigma.is_contractible() == CONTRACTIBLE;
-  return res;
+  return link_blocker_sigma.is_contractible() == CONTRACTIBLE;
 }
 
 /**
@@ -133,7 +131,7 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_all_popable_blockers(Ve
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle v) {
   // we remove the blockers that are not consistent anymore
-  update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(v));
+  update_blockers_after_remove_star_of_vertex_or_edge(Simplex(v));
   while (this->degree(v) > 0) {
     Vertex_handle w(* (adjacent_vertices(v.vertex, this->skeleton).first));
     this->remove_edge(v, w);
@@ -148,7 +146,7 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle v) {
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::update_blockers_after_remove_star_of_vertex_or_edge(
-    const Simplex_handle& simplex_to_be_removed) {
+    const Simplex& simplex_to_be_removed) {
   std::list <Blocker_handle> blockers_to_update;
   if (simplex_to_be_removed.empty()) return;
 
@@ -159,7 +157,7 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::update_blockers_after_remove_s
   }
 
   for (auto blocker_to_update : blockers_to_update) {
-    Simplex_handle sub_blocker_to_be_added;
+    Simplex sub_blocker_to_be_added;
     bool sub_blocker_need_to_be_added =
         (blocker_to_update->dimension() - simplex_to_be_removed.dimension()) >= 2;
     if (sub_blocker_need_to_be_added) {
@@ -178,7 +176,7 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::update_blockers_after_remove_s
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle a, Vertex_handle b) {
-  update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(a, b));
+  update_blockers_after_remove_star_of_vertex_or_edge(Simplex(a, b));
   // we remove the edge
   this->remove_edge(a, b);
 }
@@ -195,7 +193,7 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Edge_handle e) {
  * Remove the star of the simplex 'sigma' which needs to belong to the complex
  */
 template<typename SkeletonBlockerDS>
-void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(const Simplex_handle& sigma) {
+void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(const Simplex& sigma) {
   assert(this->contains(sigma));
   if (sigma.dimension() == 0) {
     remove_star(sigma.first_vertex());
@@ -207,34 +205,41 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(const Simplex_hand
   }
 }
 
-/**
- * @brief add a maximal simplex plus all its cofaces. All vertices lower than the higher vertex of
- * sigma must already be present.
- * @details the simplex must have dimension greater than one (otherwise use add_vertex or add_edge).
- */
 template<typename SkeletonBlockerDS>
-void Skeleton_blocker_complex<SkeletonBlockerDS>::add_simplex(const Simplex_handle& sigma) {
+void Skeleton_blocker_complex<SkeletonBlockerDS>::add_simplex(const Simplex& sigma) {
+  // to add a simplex s, all blockers included in s are first removed
+  // and then all simplex in the coboundary of s are added as blockers
   assert(!this->contains(sigma));
   assert(sigma.dimension() > 1);
+  if (!contains_vertices(sigma)) {
+    std::cerr << "add_simplex: Some vertices were not present in the complex, adding them" << std::endl;
+    size_t num_vertices_to_add = sigma.last_vertex() - this->num_vertices() + 1;
+    for (size_t i = 0; i < num_vertices_to_add; ++i)
+      this->add_vertex();
+  }
+  assert(contains_vertices(sigma));
+  if (!contains_edges(sigma))
+    add_edge(sigma);
+  remove_blocker_include_in_simplex(sigma);
+  add_blockers_after_simplex_insertion(sigma);
+}
 
-  int num_vertex_to_add = 0;
-  for (auto v : sigma)
-    if (!contains_vertex(v)) ++num_vertex_to_add;
-  while (num_vertex_to_add--) add_vertex();
 
-  for (auto u_it = sigma.begin(); u_it != sigma.end(); ++u_it)
-    for (auto v_it = u_it; ++v_it != sigma.end(); /**/) {
-      // std::cout << "add edge" << *u_it << " " << *v_it << std::endl;
-      add_edge(*u_it, *v_it);
-    }
-  remove_blocker_include_in_simplex(sigma);
+
+template<typename SkeletonBlockerDS>
+void Skeleton_blocker_complex<SkeletonBlockerDS>::add_blockers_after_simplex_insertion(Simplex sigma) {
+  if (sigma.dimension() < 1) return;
+
+  for (auto s : coboundary_range(sigma)) {
+    this->add_blocker(s);
+  }
 }
 
 /**
  * remove all blockers that contains sigma
  */
 template<typename SkeletonBlockerDS>
-void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_containing_simplex(const Simplex_handle& sigma) {
+void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_containing_simplex(const Simplex& sigma) {
   std::vector <Blocker_handle> blockers_to_remove;
   for (auto blocker : this->blocker_range(sigma.first_vertex())) {
     if (blocker->contains(sigma))
@@ -248,14 +253,26 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_containing_simp
  * remove all blockers that contains sigma
  */
 template<typename SkeletonBlockerDS>
-void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_include_in_simplex(const Simplex_handle& sigma) {
-  std::vector <Blocker_handle> blockers_to_remove;
-  for (auto blocker : this->blocker_range(sigma.first_vertex())) {
-    if (sigma.contains(*blocker))
-      blockers_to_remove.push_back(blocker);
+void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_include_in_simplex(const Simplex& sigma) {
+  // TODO(DS): write efficiently by using only superior blockers
+  // eg for all s, check blockers whose vertices are all greater than s
+  std::set <Blocker_handle> blockers_to_remove;
+  for (auto s : sigma) {
+    for (auto blocker : this->blocker_range(s)) {
+      if (sigma.contains(*blocker))
+        blockers_to_remove.insert(blocker);
+    }
   }
-  for (auto blocker_to_update : blockers_to_remove)
+  for (auto blocker_to_update : blockers_to_remove) {
+    auto s = *blocker_to_update;
     this->delete_blocker(blocker_to_update);
+    // now if there is a vertex v in the link of s
+    // and v is not included in sigma then v.s is a blocker
+    // (all faces of v.s are there since v belongs to the link of s)
+    for (const auto& b : coboundary_range(s))
+      if (!sigma.contains(b))
+        this->add_blocker(b);
+  }
 }
 
 /**
@@ -265,21 +282,21 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_include_in_simp
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::tip_blockers(Vertex_handle a, Vertex_handle b,
-                                                               std::vector<Simplex_handle> & buffer) const {
+                                                               std::vector<Simplex> & buffer) const {
   for (auto const & blocker : this->const_blocker_range(a)) {
-    Simplex_handle beta = (*blocker);
+    Simplex beta = (*blocker);
     beta.remove_vertex(a);
     buffer.push_back(beta);
   }
 
-  Simplex_handle n;
+  Simplex n;
   this->add_neighbours(b, n);
   this->remove_neighbours(a, n);
   n.remove_vertex(a);
 
 
   for (Vertex_handle y : n) {
-    Simplex_handle beta;
+    Simplex beta;
     beta.add_vertex(y);
     buffer.push_back(beta);
   }
@@ -296,7 +313,7 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::tip_blockers(Vertex_handle a,
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x) {
-  this->add_edge(a, x);
+  this->add_edge_without_blockers(a, x);
   if (this->visitor) this->visitor->on_swaped_edge(a, b, x);
   this->remove_edge(b, x);
 }
@@ -341,13 +358,13 @@ Skeleton_blocker_complex<SkeletonBlockerDS>::contract_edge(Vertex_handle a, Vert
   assert(this->contains_vertex(b));
 
   if (this->contains_edge(a, b))
-    this->add_edge(a, b);
+    this->add_edge_without_blockers(a, b);
 
   // if some blockers passes through 'ab', we need to remove them.
   if (!link_condition(a, b))
     delete_blockers_around_edge(a, b);
 
-  std::set<Simplex_handle> blockers_to_add;
+  std::set<Simplex> blockers_to_add;
 
   get_blockers_to_be_added_after_contraction(a, b, blockers_to_add);
 
@@ -367,9 +384,8 @@ Skeleton_blocker_complex<SkeletonBlockerDS>::contract_edge(Vertex_handle a, Vert
 }
 
 template<typename SkeletonBlockerDS>
-void
-Skeleton_blocker_complex<SkeletonBlockerDS>::get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, 
-    std::set<Simplex_handle>& blockers_to_add) {
+void Skeleton_blocker_complex<SkeletonBlockerDS>::get_blockers_to_be_added_after_contraction(Vertex_handle a,
+    Vertex_handle b, std::set<Simplex>& blockers_to_add) {
   blockers_to_add.clear();
 
   typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > LinkComplexType;
@@ -377,19 +393,19 @@ Skeleton_blocker_complex<SkeletonBlockerDS>::get_blockers_to_be_added_after_cont
   LinkComplexType link_a(*this, a);
   LinkComplexType link_b(*this, b);
 
-  std::vector<Simplex_handle> vector_alpha, vector_beta;
+  std::vector<Simplex> vector_alpha, vector_beta;
 
   tip_blockers(a, b, vector_alpha);
   tip_blockers(b, a, vector_beta);
 
   for (auto alpha = vector_alpha.begin(); alpha != vector_alpha.end(); ++alpha) {
     for (auto beta = vector_beta.begin(); beta != vector_beta.end(); ++beta) {
-      Simplex_handle sigma = *alpha;
+      Simplex sigma = *alpha;
       sigma.union_vertices(*beta);
       Root_simplex_handle sigma_id = this->get_id(sigma);
       if (this->contains(sigma) &&
           proper_faces_in_union<Skeleton_blocker_complex < SkeletonBlockerDS >> (sigma_id, link_a, link_b)) {
-        // Blocker_handle blocker = new Simplex_handle(sigma);
+        // Blocker_handle blocker = new Simplex(sigma);
         sigma.add_vertex(a);
         blockers_to_add.insert(sigma);
       }