skbl correction bug constructor + contraction add garland heckbert example

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@466 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 1789e903625df3c4e7c689fa4888bebd86e616eb

5 files changed, 1609 insertions, 1504 deletions

diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h
index bc832d5d..aaa682c3 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h
@@ -1,24 +1,24 @@
- /*    This file is part of the Gudhi Library. The Gudhi library
-  *    (Geometric Understanding in Higher Dimensions) is a generic C++
-  *    library for computational topology.
-  *
-  *    Author(s):       David Salinas
-  *
-  *    Copyright (C) 2014  INRIA Sophia Antipolis-Mediterranee (France)
-  *
-  *    This program is free software: you can redistribute it and/or modify
-  *    it under the terms of the GNU General Public License as published by
-  *    the Free Software Foundation, either version 3 of the License, or
-  *    (at your option) any later version.
-  *
-  *    This program is distributed in the hope that it will be useful,
-  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
-  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  *    GNU General Public License for more details.
-  *
-  *    You should have received a copy of the GNU General Public License
-  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
-  */
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       David Salinas
+ *
+ *    Copyright (C) 2014  INRIA Sophia Antipolis-Mediterranee (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
 #ifndef SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_SKELETON_BLOCKER_OFF_IO_H_
 #define SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_SKELETON_BLOCKER_OFF_IO_H_
 
@@ -35,72 +35,173 @@ namespace skbl {
  *@brief Off reader visitor that can be passed to Off_reader to read a Skeleton_blocker_complex.
  */
 template<typename Complex>
-class Skeleton_blocker_off_visitor_reader {
-  Complex& complex_;
-  typedef typename Complex::Vertex_handle Vertex_handle;
-  typedef typename Complex::Point Point;
+class Skeleton_blocker_off_flag_visitor_reader {
+	Complex& complex_;
+	typedef typename Complex::Vertex_handle Vertex_handle;
+	typedef typename Complex::Point Point;
 
-  const bool load_only_points_;
+	const bool load_only_points_;
 
- public:
-  explicit Skeleton_blocker_off_visitor_reader(Complex& complex, bool load_only_points = false):
-    complex_(complex),
-    load_only_points_(load_only_points) {}
+public:
+	explicit Skeleton_blocker_off_flag_visitor_reader(Complex& complex, bool load_only_points = false):
+	complex_(complex),
+	load_only_points_(load_only_points) {}
 
 
-  void init(int dim, int num_vertices, int num_faces, int num_edges) {
-    // todo do an assert to check that this number are correctly read
-    // todo reserve size for vector points
-  }
+	void init(int dim, int num_vertices, int num_faces, int num_edges) {
+		// todo do an assert to check that this number are correctly read
+		// todo reserve size for vector points
+	}
 
 
-  void point(const std::vector<double>& point) {
-    complex_.add_vertex(point);
-  }
+	void point(const std::vector<double>& point) {
+		complex_.add_vertex(Point(point.begin(),point.end()));
+	}
 
-  void maximal_face(const std::vector<int>& face) {
-    if (!load_only_points_) {
-      for (size_t i = 0; i < face.size(); ++i)
-        for (size_t j = i+1; j < face.size(); ++j) {
-          complex_.add_edge(Vertex_handle(face[i]), Vertex_handle(face[j]));
-        }
-      }
-  }
+	void maximal_face(const std::vector<int>& face) {
+		if (!load_only_points_) {
+			for (size_t i = 0; i < face.size(); ++i)
+				for (size_t j = i+1; j < face.size(); ++j) {
+					complex_.add_edge(Vertex_handle(face[i]), Vertex_handle(face[j]));
+				}
+		}
+	}
 
-  void done() {}
+	void done() {}
 };
 
 /**
-*@brief Class that allows to load a Skeleton_blocker_complex from an off file.
-*/
+ *@brief Off reader visitor that can be passed to Off_reader to read a Skeleton_blocker_complex.
+ */
+template<typename Complex>
+class Skeleton_blocker_off_visitor_reader {
+	Complex& complex_;
+	typedef typename Complex::Vertex_handle Vertex_handle;
+	typedef typename Complex::Simplex_handle Simplex_handle;
+	typedef typename Complex::Point Point;
+
+	const bool load_only_points_;
+	std::vector<Point> points_;
+	std::vector<Simplex_handle> maximal_faces_;
+
+public:
+	explicit Skeleton_blocker_off_visitor_reader(Complex& complex, bool load_only_points = false):
+	complex_(complex),
+	load_only_points_(load_only_points) {}
+
+
+	void init(int dim, int num_vertices, int num_faces, int num_edges){
+		maximal_faces_.reserve(num_faces);
+		points_.reserve(num_vertices);
+	}
+
+
+	void point(const std::vector<double>& point) {
+		points_.emplace_back(point.begin(),point.end());
+	}
+
+	void maximal_face(const std::vector<int>& face) {
+		if (!load_only_points_) {
+			Simplex_handle s;
+			for(auto x : face)
+				s.add_vertex(Vertex_handle(x));
+			maximal_faces_.emplace_back(s);
+		}
+	}
+
+	void done() {
+		complex_ = make_complex_from_top_faces(
+				maximal_faces_.begin(),
+				maximal_faces_.end(),
+				points_.begin(),
+				points_.end()
+		);
+	}
+};
+
+
+/**
+ *@brief Class that allows to load a Skeleton_blocker_complex from an off file.
+ */
 template<typename Complex>
 class Skeleton_blocker_off_reader {
- public:
-  /**
-   * name_file : file to read
-   * read_complex : complex that will receive the file content
-   * read_only_points : specify true if only the points must be read
-   */
-  Skeleton_blocker_off_reader(const std::string & name_file, Complex& read_complex, bool read_only_points = false):valid_(false) {
-    std::ifstream stream(name_file);
-    if (stream.is_open()) {
-      Skeleton_blocker_off_visitor_reader<Complex> off_visitor(read_complex, read_only_points);
-      Off_reader off_reader(stream);
-      valid_ = off_reader.read(off_visitor);
-    }
-  }
-
-  /**
-   * return true iff reading did not meet problems.
-   */
-  bool is_valid() const {
-    return valid_;
-  }
-
- private:
-  bool valid_;
+public:
+	/**
+	 * name_file : file to read
+	 * read_complex : complex that will receive the file content
+	 * read_only_points : specify true if only the points must be read
+	 */
+	Skeleton_blocker_off_reader(const std::string & name_file, Complex& read_complex, bool read_only_points = false,bool is_flag = false):valid_(false) {
+		std::ifstream stream(name_file);
+		if (stream.is_open()) {
+			if(is_flag){
+				Skeleton_blocker_off_flag_visitor_reader<Complex> off_visitor(read_complex, read_only_points);
+				Off_reader off_reader(stream);
+				valid_ = off_reader.read(off_visitor);
+			}
+			else{
+				Skeleton_blocker_off_visitor_reader<Complex> off_visitor(read_complex, read_only_points);
+				Off_reader off_reader(stream);
+				valid_ = off_reader.read(off_visitor);
+			}
+		}
+	}
+
+	/**
+	 * return true iff reading did not meet problems.
+	 */
+	bool is_valid() const {
+		return valid_;
+	}
+
+private:
+	bool valid_;
+};
+
+
+template<typename Complex>
+class Skeleton_blocker_off_writer{
+public:
+	/**
+	 * name_file : file where the off will be written
+	 * save_complex : complex that be outputted in the file
+	 * for now only save triangles.
+	 */
+	Skeleton_blocker_off_writer(const std::string & name_file, const Complex& save_complex){
+		typedef typename Complex::Vertex_handle Vertex_handle;
+
+		std::ofstream stream(name_file);
+		if (stream.is_open()) {
+			stream<<"OFF\n";
+			size_t num_triangles = std::distance(save_complex.triangle_range().begin(),save_complex.triangle_range().end());
+			stream<<save_complex.num_vertices()<<" "<<num_triangles<< " 0 \n";
+
+			//in case the complex has deactivated some vertices, eg only has vertices 0 2 5 7 for instance
+			//we compute a map from 0 2 5 7 to 0 1 2 3
+			std::map<Vertex_handle,size_t> vertex_num;
+			size_t current_vertex=0;
+
+			for(auto v : save_complex.vertex_range()){
+				vertex_num[v]=current_vertex++;
+				const auto& pt(save_complex.point(v));
+				for(auto it = pt.begin();it!=pt.end();++it)
+					stream<<*it<<" ";
+				stream<<std::endl;
+			}
+
+			for(const auto & t : save_complex.triangle_range()){
+				stream<<"3 ";
+				for(auto x : t)
+					stream<<vertex_num[x]<<" ";
+				stream<<std::endl;
+			}
+			stream.close();
+		} else std::cerr<<"could not open file "<<name_file<<std::endl;
+	}
+
 };
 
+
 }  // namespace skbl
 
 
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h
index 717b6052..434538fe 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h
@@ -64,9 +64,6 @@ class Skeleton_blocker_simplex {
   //@{
 
   // Skeleton_blocker_simplex():simplex_set() {}
-  /**
-   * Clear the simplex
-   */
   void clear() {
     simplex_set.clear();
   }
@@ -121,7 +118,6 @@ class Skeleton_blocker_simplex {
    * Add the vertex v to the simplex:
    * Note that adding two times the same vertex is
    * the same that adding it once.
-   * \f$ (*this) \leftarrow (*this) \cup \{ v \} \f$
    */
   void add_vertex(T v) {
     simplex_set.insert(v);
@@ -129,15 +125,13 @@ class Skeleton_blocker_simplex {
 
   /**
    * Remove the vertex v from the simplex:
-   * \f$ (*this) \leftarrow (*this) \setminus \{ v \} \f$
    */
   void remove_vertex(T v) {
     simplex_set.erase(v);
   }
 
   /**
-   * Intersects the simplex with the simplex a:
-   * \f$ (*this) \leftarrow (*this) \cap a \f$
+   * Intersects the simplex with the simplex a.
    */
   void intersection(const Skeleton_blocker_simplex & a) {
     std::vector<T> v;
@@ -152,8 +146,7 @@ class Skeleton_blocker_simplex {
   }
 
   /**
-   * Substracts a from the simplex:
-   * \f$ (*this) \leftarrow (*this) \setminus a \f$
+   * Substracts a from the simplex.
    */
   void difference(const Skeleton_blocker_simplex & a) {
     std::vector<T> v;
@@ -169,8 +162,7 @@ class Skeleton_blocker_simplex {
   }
 
   /**
-   * Add vertices of a to the simplex:
-   * \f$ (*this) \leftarrow (*this) \cup a \f$
+   * Add vertices of a to the simplex.
    */
   void union_vertices(const Skeleton_blocker_simplex & a) {
     std::vector<T> v;
@@ -234,7 +226,7 @@ class Skeleton_blocker_simplex {
     return *(simplex_set.rbegin());
   }
   /**
-   * @return true iff the simplex contains the simplex a, i.e. iff \f$ a \subset (*this) \f$.
+   * @return true iff the simplex contains the simplex a.
    */
   bool contains(const Skeleton_blocker_simplex & a) const {
     return includes(simplex_set.cbegin(), simplex_set.cend(),
@@ -278,7 +270,7 @@ class Skeleton_blocker_simplex {
     auto last2 = a.end();
 
     while (first2 != last2) {
-      // we ignore vertices of b
+      // we ignore vertices x
       if (x == *first2) {
         ++first2;
       } else {
@@ -293,7 +285,7 @@ class Skeleton_blocker_simplex {
   }
 
   /**
-   * @return true iff the simplex contains the difference \f$ a \setminus \{ x \} \f$.
+   * @return true iff the simplex contains the difference \f$ a \setminus \{ x,y \} \f$.
    */
   bool contains_difference(const Skeleton_blocker_simplex& a, T x, T y) const {
     auto first1 = begin();
@@ -303,7 +295,7 @@ class Skeleton_blocker_simplex {
     auto last2 = a.end();
 
     while (first2 != last2) {
-      // we ignore vertices of b
+      // we ignore vertices of x,y
       if (x == *first2 || y == *first2) {
         ++first2;
       } else {
@@ -317,16 +309,10 @@ class Skeleton_blocker_simplex {
     return true;
   }
 
-  /**
-   * @return true iff the simplex contains the vertex v, i.e. iff \f$ v \in (*this) \f$.
-   */
   bool contains(T v) const {
     return (simplex_set.find(v) != simplex_set.end());
   }
 
-  /**
-   * @return \f$ (*this) \cap a = \emptyset \f$.
-   */
   bool disjoint(const Skeleton_blocker_simplex& a) const {
     std::vector<T> v;
     v.reserve(std::min(simplex_set.size(), a.simplex_set.size()));
@@ -354,9 +340,6 @@ class Skeleton_blocker_simplex {
 
   //@}
 
-  /**
-   * Display a simplex
-   */
   friend std::ostream& operator <<(std::ostream& o,
                                    const Skeleton_blocker_simplex & sigma) {
     bool first = true;
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
index 7b7ac7a9..81fc28cb 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
@@ -63,977 +63,988 @@ namespace skbl {
  */
 template<class SkeletonBlockerDS>
 class Skeleton_blocker_complex {
-  template<class ComplexType> friend class Complex_vertex_iterator;
-  template<class ComplexType> friend class Complex_neighbors_vertices_iterator;
-  template<class ComplexType> friend class Complex_edge_iterator;
-  template<class ComplexType> friend class Complex_edge_around_vertex_iterator;
-
-  template<class ComplexType> friend class Skeleton_blocker_link_complex;
-  template<class ComplexType> friend class Skeleton_blocker_link_superior;
-  template<class ComplexType> friend class Skeleton_blocker_sub_complex;
-
- public:
-  /**
-   * @brief The type of stored vertex node, specified by the template SkeletonBlockerDS
-   */
-  typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex;
-
-  /**
-   * @brief The type of stored  edge node, specified by the template SkeletonBlockerDS
-   */
-  typedef typename SkeletonBlockerDS::Graph_edge Graph_edge;
-
-  typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle;
-
-  /**
-   * @brief The type of an handle to a vertex of the complex.
-   */
-  typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle;
-  typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle;
-
-  /**
-   * @brief A ordered set of integers that represents a simplex.
-   */
-  typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle;
-  typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle;
-
-  /**
-   * @brief Handle to a blocker of the complex.
-   */
-  typedef Simplex_handle* Blocker_handle;
-
-  typedef typename Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_iterator;
-  typedef typename Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator;
-
- protected:
-  typedef typename boost::adjacency_list<boost::setS,  // edges
-  boost::vecS,  // vertices
-  boost::undirectedS, Graph_vertex, Graph_edge> Graph;
-  // todo/remark : edges are not sorted, it heavily penalizes computation for SuperiorLink
-  // (eg Link with greater vertices)
-  // that burdens simplex iteration / complex initialization via list of simplices.
-  // to avoid that, one should modify the graph by storing two lists of adjacency for every
-  // vertex, the one with superior and the one with lower vertices, that way there is
-  // no more extra cost for computation of SuperiorLink
-  typedef typename boost::graph_traits<Graph>::vertex_iterator boost_vertex_iterator;
-  typedef typename boost::graph_traits<Graph>::edge_iterator boost_edge_iterator;
-
- protected:
-  typedef typename boost::graph_traits<Graph>::adjacency_iterator boost_adjacency_iterator;
-
- public:
-  /**
-   * @brief Handle to an edge of the complex.
-   */
-  typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle;
-
- protected:
-  typedef std::multimap<Vertex_handle, Simplex_handle *> BlockerMap;
-  typedef typename std::multimap<Vertex_handle, Simplex_handle *>::value_type BlockerPair;
-  typedef typename std::multimap<Vertex_handle, Simplex_handle *>::iterator BlockerMapIterator;
-  typedef typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator BlockerMapConstIterator;
-
- protected:
-  int num_vertices_;
-  int num_blockers_;
-
-  typedef Skeleton_blocker_complex_visitor<Vertex_handle> Visitor;
-  // typedef Visitor* Visitor_ptr;
-  Visitor* visitor;
-
-  /**
-   * @details If 'x' is a Vertex_handle of a vertex in the complex then degree[x] = d is its degree.
-   *
-   * This quantity is updated when adding/removing edge.
-   *
-   * This is useful because the operation
-   * list.size() is done in linear time.
-   */
-  std::vector<boost_vertex_handle> degree_;
-  Graph skeleton; /** 1-skeleton of the simplicial complex. */
-
-  /** Each vertex can access to the blockers passing through it. */
-  BlockerMap blocker_map_;
-
- public:
-  /////////////////////////////////////////////////////////////////////////////
-  /** @name Constructors, Destructors
-   */
-  //@{
-
-  /**
-   *@brief constructs a simplicial complex with a given number of vertices and a visitor.
-   */
-  explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL)
-      : visitor(visitor_) {
-    clear();
-    for (int i = 0; i < num_vertices_; ++i) {
-      add_vertex();
-    }
-  }
-
- private:
-  typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie;
-
-  template<typename SimpleHandleOutputIterator>
-  /**
-   * return a vector of simplex trie for every vertex
-   */
-  std::vector<STrie*> compute_cofaces(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) {
-    std::vector<STrie*> cofaces(num_vertices(), 0);
-    for (auto i = 0u; i < num_vertices(); ++i)
-      cofaces[i] = new STrie(Vertex_handle(i));
-    for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) {
-      if (s_it->dimension() >= 1)
-        cofaces[s_it->first_vertex()]->add_simplex(*s_it);
-    }
-    return cofaces;
-  }
-
-
- public:
-  /**
-   * @brief Constructor with a list of simplices.
-   * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not).
-   */
-  template<typename SimpleHandleOutputIterator>
-  Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, bool is_flag_complex = false, Visitor* visitor_ = NULL)
-      : num_vertices_(0),
-      num_blockers_(0),
-      visitor(visitor_) {
-    std::vector<std::pair<Vertex_handle, Vertex_handle>> edges;
-
-    // first pass to add vertices and edges
-    int num_vertex = -1;
-    for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) {
-      if (s_it->dimension() == 0) num_vertex = (std::max)(num_vertex, s_it->first_vertex().vertex);
-      if (s_it->dimension() == 1) edges.emplace_back(s_it->first_vertex(), s_it->last_vertex());
-    }
-    while (num_vertex-- >= 0) add_vertex();
-
-    for (const auto& e : edges)
-      add_edge(e.first, e.second);
-
-    if (!is_flag_complex) {
-      // need to compute blockers
-      std::vector<STrie*> cofaces(compute_cofaces(simplex_begin, simplex_end));
-      std::vector<Simplex_handle> max_faces;
-
-      for (auto i = 0u; i < num_vertices(); ++i) {
-        auto max_faces_i = cofaces[i]->maximal_faces();
-        max_faces.insert(max_faces.end(), max_faces_i.begin(), max_faces_i.end());
-      }
-
-      // all time here
-
-      // for every maximal face s, it picks its first vertex v and check for all nv \in N(v)
-      // if s union nv is in the complex, if not it is a blocker.
-      for (auto max_face : max_faces) {
-        if (max_face.dimension() > 0) {
-          auto first_v = max_face.first_vertex();
-          for (auto nv : vertex_range(first_v)) {
-            if (!max_face.contains(nv) && max_face.first_vertex() < nv) {
-              // check that all edges in vertices(max_face)\cup nv are here
-              // since max_face is a simplex, we only need to check that edges
-              // (x nv) where x \in max_face are present
-              bool all_edges_here = true;
-              for (auto x : max_face)
-                if (!contains_edge(x, nv)) {
-                  all_edges_here = false;
-                  break;
-                }
-              if (all_edges_here) {  // eg  this->contains(max_face)
-                max_face.add_vertex(nv);
-                if (!cofaces[first_v]->contains(max_face)) {
-                  // if there exists a blocker included in max_face, we remove it
-                  // as it is not a minimum missing face
-                  // the other alternative would be to check to all properfaces
-                  // are in the complex before adding a blocker but that
-                  // would be more expensive if there are few blockers
-                  std::vector<Blocker_handle> blockers_to_remove;
-                  for (auto b : blocker_range(first_v))
-                    if (b->contains(max_face))
-                      blockers_to_remove.push_back(b);
-                  for (auto b : blockers_to_remove)
-                    this->delete_blocker(b);
-
-                  add_blocker(max_face);
-                }
-                max_face.remove_vertex(nv);
-              }
-            }
-          }
-        }
-      }
-      for (auto i = 0u; i < num_vertices(); ++i)
-        delete(cofaces[i]);
-    }
-  }
-
-  // We cannot use the default copy constructor since we need
-  // to make a copy of each of the blockers
-
-  Skeleton_blocker_complex(const Skeleton_blocker_complex& copy) {
-    visitor = NULL;
-    degree_ = copy.degree_;
-    skeleton = Graph(copy.skeleton);
-    num_vertices_ = copy.num_vertices_;
-
-    num_blockers_ = 0;
-    // we copy the blockers
-    for (auto blocker : copy.const_blocker_range()) {
-      add_blocker(*blocker);
-    }
-  }
-
-  /**
-   */
-  Skeleton_blocker_complex& operator=(const Skeleton_blocker_complex& copy) {
-    clear();
-    visitor = NULL;
-    degree_ = copy.degree_;
-    skeleton = Graph(copy.skeleton);
-    num_vertices_ = copy.num_vertices_;
-
-    num_blockers_ = 0;
-    // we copy the blockers
-    for (auto blocker : copy.const_blocker_range())
-      add_blocker(*blocker);
-    return *this;
-  }
-
-
-  /**
-   * return true if both complexes have the same simplices.
-   */
-  bool operator==(const Skeleton_blocker_complex& other) const{
-	  if(other.num_vertices() != num_vertices()) return false;
-	  if(other.num_edges() != num_edges()) return false;
-	  if(other.num_blockers() != num_blockers()) return false;
-
-	  for(auto v : vertex_range())
-		  if(!other.contains_vertex(v)) return false;
-
-	  for(auto e : edge_range())
-		  if(!other.contains_edge(first_vertex(e),second_vertex(e))) return false;
-
-	  for(const auto b : const_blocker_range())
-		  if(!other.contains_blocker(*b)) return false;
-
-	  return true;
-  }
-
-  bool operator!=(const Skeleton_blocker_complex& other) const{
-	  return !(*this == other);
-  }
-
-  /**
-   * The destructor delete all blockers allocated.
-   */
-  virtual ~Skeleton_blocker_complex() {
-    clear();
-  }
-
-  /**
-   * @details Clears the simplicial complex. After a call to this function,
-   * blockers are destroyed. The 1-skeleton and the set of blockers
-   * are both empty.
-   */
-  virtual void clear() {
-    // xxx for now the responsabilty of freeing the visitor is for
-    // the user
-    visitor = NULL;
-
-    degree_.clear();
-    num_vertices_ = 0;
-
-    remove_blockers();
-
-    skeleton.clear();
-  }
-
-  /**
-   *@brief allows to change the visitor.
-   */
-  void set_visitor(Visitor* other_visitor) {
-    visitor = other_visitor;
-  }
-
-  //@}
-
-  /////////////////////////////////////////////////////////////////////////////
-  /** @name Vertices operations
-   */
-  //@{
- public:
-  /**
-   * @brief Return a local Vertex_handle of a vertex given a global one.
-   * @remark Assume that the vertex is present in the complex.
-   */
-  Vertex_handle operator[](Root_vertex_handle global) const {
-    auto local(get_address(global));
-    assert(local);
-    return *local;
-  }
-
-  /**
-   * @brief Return the vertex node associated to local Vertex_handle.
-   * @remark Assume that the vertex is present in the complex.
-   */
-  Graph_vertex& operator[](Vertex_handle address) {
-    assert(
-           0 <= address.vertex && address.vertex < boost::num_vertices(skeleton));
-    return skeleton[address.vertex];
-  }
-
-  /**
-   * @brief Return the vertex node associated to local Vertex_handle.
-   * @remark Assume that the vertex is present in the complex.
-   */
-  const Graph_vertex& operator[](Vertex_handle address) const {
-    assert(
-           0 <= address.vertex && address.vertex < boost::num_vertices(skeleton));
-    return skeleton[address.vertex];
-  }
-
-  /**
-   * @brief Adds a vertex to the simplicial complex and returns its Vertex_handle.
-   */
-  Vertex_handle add_vertex() {
-    Vertex_handle address(boost::add_vertex(skeleton));
-    num_vertices_++;
-    (*this)[address].activate();
-    // safe since we now that we are in the root complex and the field 'address' and 'id'
-    // are identical for every vertices
-    (*this)[address].set_id(Root_vertex_handle(address.vertex));
-    degree_.push_back(0);
-    if (visitor)
-      visitor->on_add_vertex(address);
-    return address;
-  }
-
-  /**
-   * @brief Remove a vertex from the simplicial complex
-   * @remark It just deactivates the vertex with a boolean flag but does not
-   * remove it from vertices from complexity issues.
-   */
-  void remove_vertex(Vertex_handle address) {
-    assert(contains_vertex(address));
-    // We remove b
-    boost::clear_vertex(address.vertex, skeleton);
-    (*this)[address].deactivate();
-    num_vertices_--;
-    degree_[address.vertex] = -1;
-    if (visitor)
-      visitor->on_remove_vertex(address);
-  }
-
-  /**
-   */
-  bool contains_vertex(Vertex_handle u) const {
-    if (u.vertex < 0 || u.vertex >= boost::num_vertices(skeleton))
-      return false;
-    return (*this)[u].is_active();
-  }
-
-  /**
-   */
-  bool contains_vertex(Root_vertex_handle u) const {
-    boost::optional<Vertex_handle> address = get_address(u);
-    return address && (*this)[*address].is_active();
-  }
-
-  /**
-   * @return true iff the simplicial complex contains all vertices
-   * of simplex sigma
-   */
-  bool contains_vertices(const Simplex_handle & sigma) const {
-    for (auto vertex : sigma)
-      if (!contains_vertex(vertex))
-        return false;
-    return true;
-  }
-
-  /**
-   * @brief Given an Id return the address of the vertex having this Id in the complex.
-   * @remark For a simplicial complex, the address is the id but it may not be the case for a SubComplex.
-   */
-  virtual boost::optional<Vertex_handle> get_address(
-                                                     Root_vertex_handle id) const {
-    boost::optional<Vertex_handle> res;
-    if (id.vertex < boost::num_vertices(skeleton))
-      res = Vertex_handle(id.vertex);  // xxx
-    return res;
-  }
-
-  /**
-   * return the id of a vertex of adress local present in the graph
-   */
-  Root_vertex_handle get_id(Vertex_handle local) const {
-    assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton));
-    return (*this)[local].get_id();
-  }
-
-  /**
-   * @brief Convert an address of a vertex of a complex to the address in
-   * the current complex.
-   * @details
-   * If the current complex is a sub (or sup) complex of 'other', it converts
-   * the address of a vertex v expressed in 'other' to the address of the vertex
-   * v in the current one.
-   * @remark this methods uses Root_vertex_handle to identify the vertex and
-   * assumes the vertex is present in the current complex.
-   */
-  Vertex_handle convert_handle_from_another_complex(
-                                                    const Skeleton_blocker_complex& other, Vertex_handle vh_in_other) const {
-    auto vh_in_current_complex = get_address(other.get_id(vh_in_other));
-    assert(vh_in_current_complex);
-    return *vh_in_current_complex;
-  }
-
-  /**
-   * @brief return the graph degree of a vertex.
-   */
-  int degree(Vertex_handle local) const {
-    assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton));
-    return degree_[local.vertex];
-  }
-
-  //@}
-
-  /////////////////////////////////////////////////////////////////////////////
-  /** @name Edges operations
-   */
-  //@{
- public:
-  /**
-   * @brief return an edge handle if the two vertices forms
-   * an edge in the complex
-   */
-  boost::optional<Edge_handle> operator[](
-                                          const std::pair<Vertex_handle, Vertex_handle>& ab) const {
-    boost::optional<Edge_handle> res;
-    std::pair<Edge_handle, bool> edge_pair(
-                                           boost::edge(ab.first.vertex, ab.second.vertex, skeleton));
-    if (edge_pair.second)
-      res = edge_pair.first;
-    return res;
-  }
-
-  /**
-   * @brief returns the stored node associated to an edge
-   */
-  Graph_edge& operator[](Edge_handle edge_handle) {
-    return skeleton[edge_handle];
-  }
-
-  /**
-   * @brief returns the stored node associated to an edge
-   */
-  const Graph_edge& operator[](Edge_handle edge_handle) const {
-    return skeleton[edge_handle];
-  }
-
-  /**
-   * @brief returns the first vertex of an edge
-   * @details it assumes that the edge is present in the complex
-   */
-  Vertex_handle first_vertex(Edge_handle edge_handle) const {
-    return static_cast<Vertex_handle> (source(edge_handle, skeleton));
-  }
-
-  /**
-   * @brief returns the first vertex of an edge
-   * @details it assumes that the edge is present in the complex
-   */
-  Vertex_handle second_vertex(Edge_handle edge_handle) const {
-    return static_cast<Vertex_handle> (target(edge_handle, skeleton));
-  }
-
-  /**
-   * @brief returns the simplex made with the two vertices of the edge
-   * @details it assumes that the edge is present in the complex
-
-   */
-  Simplex_handle get_vertices(Edge_handle edge_handle) const {
-    auto edge((*this)[edge_handle]);
-    return Simplex_handle((*this)[edge.first()], (*this)[edge.second()]);
-  }
-
-  /**
-   * @brief Adds an edge between vertices a and b and all its cofaces.
-   */
-  Edge_handle add_edge(Vertex_handle a, Vertex_handle b) {
-    assert(contains_vertex(a) && contains_vertex(b));
-    assert(a != b);
-
-    auto edge_handle((*this)[std::make_pair(a, b)]);
-    // std::pair<Edge_handle,bool> pair_descr_bool = (*this)[std::make_pair(a,b)];
-    // Edge_handle edge_descr;
-    // bool edge_present = pair_descr_bool.second;
-    if (!edge_handle) {
-      edge_handle = boost::add_edge(a.vertex, b.vertex, skeleton).first;
-      (*this)[*edge_handle].setId(get_id(a), get_id(b));
-      degree_[a.vertex]++;
-      degree_[b.vertex]++;
-      if (visitor)
-        visitor->on_add_edge(a, b);
-    }
-    return *edge_handle;
-  }
-
-  /**
-   * @brief Adds all edges and their cofaces of a simplex to the simplicial complex.
-   */
-  void add_edges(const Simplex_handle & sigma) {
-    Simplex_handle_iterator i, j;
-    for (i = sigma.begin(); i != sigma.end(); ++i)
-      for (j = i, j++; j != sigma.end(); ++j)
-        add_edge(*i, *j);
-  }
-
-  /**
-   * @brief Removes an edge from the simplicial complex and all its cofaces.
-   * @details returns the former Edge_handle representing the edge
-   */
-  virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b) {
-    bool found;
-    Edge_handle edge;
-    tie(edge, found) = boost::edge(a.vertex, b.vertex, skeleton);
-    if (found) {
-      if (visitor)
-        visitor->on_remove_edge(a, b);
-      // if (heapCollapse.Contains(edge)) heapCollapse.Delete(edge);
-      boost::remove_edge(a.vertex, b.vertex, skeleton);
-      degree_[a.vertex]--;
-      degree_[b.vertex]--;
-    }
-    return edge;
-  }
-
-  /**
-   * @brief Removes edge and its cofaces from the simplicial complex.
-   */
-  void remove_edge(Edge_handle edge) {
-    assert(contains_vertex(first_vertex(edge)));
-    assert(contains_vertex(second_vertex(edge)));
-    remove_edge(first_vertex(edge), second_vertex(edge));
-  }
-
-  /**
-   * @brief The complex is reduced to its set of vertices.
-   * All the edges and blockers are removed.
-   */
-  void keep_only_vertices() {
-    remove_blockers();
-
-    for (auto u : vertex_range()) {
-      while (this->degree(u) > 0) {
-        Vertex_handle v(*(adjacent_vertices(u.vertex, this->skeleton).first));
-        this->remove_edge(u, v);
-      }
-    }
-  }
-
-  /**
-   * @return true iff the simplicial complex contains an edge between
-   * vertices a and b
-   */
-  bool contains_edge(Vertex_handle a, Vertex_handle b) const {
-    // if (a.vertex<0 || b.vertex <0) return false;
-    return boost::edge(a.vertex, b.vertex, skeleton).second;
-  }
-
-  /**
-   * @return true iff the simplicial complex contains all vertices
-   * and all edges of simplex sigma
-   */
-  bool contains_edges(const Simplex_handle & sigma) const {
-    for (auto i = sigma.begin(); i != sigma.end(); ++i) {
-      if (!contains_vertex(*i))
-        return false;
-      for (auto j = i; ++j != sigma.end();) {
-        if (!contains_edge(*i, *j))
-          return false;
-      }
-    }
-    return true;
-  }
-  //@}
-
-  /////////////////////////////////////////////////////////////////////////////
-  /** @name Blockers operations
-   */
-  //@{
-
-  /**
-   * @brief Adds the simplex to the set of blockers and
-   * returns a Blocker_handle toward it if was not present before and 0 otherwise.
-   */
-  Blocker_handle add_blocker(const Simplex_handle& blocker) {
-    assert(blocker.dimension() > 1);
-    if (contains_blocker(blocker)) {
-      // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl;
-      return 0;
-    } else {
-      if (visitor)
-        visitor->on_add_blocker(blocker);
-      Blocker_handle blocker_pt = new Simplex_handle(blocker);
-      num_blockers_++;
-      auto vertex = blocker_pt->begin();
-      while (vertex != blocker_pt->end()) {
-        blocker_map_.insert(BlockerPair(*vertex, blocker_pt));
-        ++vertex;
-      }
-      return blocker_pt;
-    }
-  }
-
- protected:
-  /**
-   * @brief Adds the simplex to the set of blockers
-   */
-  void add_blocker(Blocker_handle blocker) {
-    if (contains_blocker(*blocker)) {
-      // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl;
-      return;
-    } else {
-      if (visitor)
-        visitor->on_add_blocker(*blocker);
-      num_blockers_++;
-      auto vertex = blocker->begin();
-      while (vertex != blocker->end()) {
-        blocker_map_.insert(BlockerPair(*vertex, blocker));
-        ++vertex;
-      }
-    }
-  }
-
- protected:
-  /**
-   * Removes sigma from the blocker map of vertex v
-   */
-  void remove_blocker(const Blocker_handle sigma, Vertex_handle v) {
-    Complex_blocker_around_vertex_iterator blocker;
-    for (blocker = blocker_range(v).begin(); blocker != blocker_range(v).end();
-         ++blocker) {
-      if (*blocker == sigma)
-        break;
-    }
-    if (*blocker != sigma) {
-      std::cerr
-          << "bug ((*blocker).second == sigma) ie try to remove a blocker not present\n";
-      assert(false);
-    } else {
-      blocker_map_.erase(blocker.current_position());
-    }
-  }
-
- public:
-  /**
-   * @brief Removes the simplex from the set of blockers.
-   * @remark sigma has to belongs to the set of blockers
-   */
-  void remove_blocker(const Blocker_handle sigma) {
-    for (auto vertex : *sigma)
-      remove_blocker(sigma, vertex);
-    num_blockers_--;
-  }
-
-  /**
-   * @brief Remove all blockers, in other words, it expand the simplicial
-   * complex to the smallest flag complex that contains it.
-   */
-  void remove_blockers() {
-    // Desallocate the blockers
-    while (!blocker_map_.empty()) {
-      delete_blocker(blocker_map_.begin()->second);
-    }
-    num_blockers_ = 0;
-    blocker_map_.clear();
-  }
-
- protected:
-  /**
-   * Removes the simplex sigma from the set of blockers.
-   * sigma has to belongs to the set of blockers
-   *
-   * @remark contrarily to delete_blockers does not call the destructor
-   */
-  void remove_blocker(const Simplex_handle& sigma) {
-    assert(contains_blocker(sigma));
-    for (auto vertex : sigma)
-      remove_blocker(sigma, vertex);
-    num_blockers_--;
-  }
-
- public:
-  /**
-   * Removes the simplex s from the set of blockers
-   * and desallocate s.
-   */
-  void delete_blocker(Blocker_handle sigma) {
-    if (visitor)
-      visitor->on_delete_blocker(sigma);
-    remove_blocker(sigma);
-    delete sigma;
-  }
-
-  /**
-   * @return true iff s is a blocker of the simplicial complex
-   */
-  bool contains_blocker(const Blocker_handle s) const {
-    if (s->dimension() < 2)
-      return false;
-
-    Vertex_handle a = s->first_vertex();
-
-    for (const auto blocker : const_blocker_range(a)) {
-      if (s == *blocker)
-        return true;
-    }
-    return false;
-  }
-
-  /**
-   * @return true iff s is a blocker of the simplicial complex
-   */
-  bool contains_blocker(const Simplex_handle & s) const {
-    if (s.dimension() < 2)
-      return false;
-
-    Vertex_handle a = s.first_vertex();
-
-    for (auto blocker : const_blocker_range(a)) {
-      if (s == *blocker)
-        return true;
-    }
-    return false;
-  }
-
- private:
-  /**
-   * @return true iff a blocker of the simplicial complex
-   * is a face of sigma.
-   */
-  bool blocks(const Simplex_handle & sigma) const {
-    for (auto blocker : const_blocker_range()) {
-      if (sigma.contains(*blocker))
-        return true;
-    }
-    return false;
-  }
-
-  //@}
-
- protected:
-  /**
-   * @details Adds to simplex the neighbours of v e.g. \f$ n \leftarrow n \cup N(v) \f$.
-   * If keep_only_superior is true then only vertices that are greater than v are added.
-   */
-  virtual void add_neighbours(Vertex_handle v, Simplex_handle & n,
-                              bool keep_only_superior = false) const {
-    boost_adjacency_iterator ai, ai_end;
-    for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end;
-         ++ai) {
-      if (keep_only_superior) {
-        if (*ai > v.vertex) {
-          n.add_vertex(Vertex_handle(*ai));
-        }
-      } else {
-        n.add_vertex(Vertex_handle(*ai));
-      }
-    }
-  }
-
-  /**
-   * @details Add to simplex res all vertices which are
-   * neighbours of alpha: ie \f$ res \leftarrow res \cup N(alpha) \f$.
-   *
-   * If 'keep_only_superior' is true then only vertices that are greater than alpha are added.
-   * todo revoir
-   *
-   */
-  virtual void add_neighbours(const Simplex_handle &alpha, Simplex_handle & res,
-                              bool keep_only_superior = false) const {
-    res.clear();
-    auto alpha_vertex = alpha.begin();
-    add_neighbours(*alpha_vertex, res, keep_only_superior);
-    for (alpha_vertex = (alpha.begin())++; alpha_vertex != alpha.end();
-         ++alpha_vertex)
-      keep_neighbours(*alpha_vertex, res, keep_only_superior);
-  }
-
-  /**
-   * @details Remove from simplex n all vertices which are
-   * not neighbours of v e.g. \f$ res \leftarrow res \cap N(v) \f$.
-   * If 'keep_only_superior' is true then only vertices that are greater than v are keeped.
-   */
-  virtual void keep_neighbours(Vertex_handle v, Simplex_handle& res,
-                               bool keep_only_superior = false) const {
-    Simplex_handle nv;
-    add_neighbours(v, nv, keep_only_superior);
-    res.intersection(nv);
-  }
-
-  /**
-   * @details Remove from simplex all vertices which are
-   * neighbours of v eg \f$ res \leftarrow res \setminus N(v) \f$.
-   * If 'keep_only_superior' is true then only vertices that are greater than v are added.
-   */
-  virtual void remove_neighbours(Vertex_handle v, Simplex_handle & res,
-                                 bool keep_only_superior = false) const {
-    Simplex_handle nv;
-    add_neighbours(v, nv, keep_only_superior);
-    res.difference(nv);
-  }
-
- public:
-  typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex> Link_complex;
-
-  /**
-   * Constructs the link of 'simplex' with points coordinates.
-   */
-  Link_complex link(Vertex_handle v) const {
-    return Link_complex(*this, Simplex_handle(v));
-  }
-
-  /**
-   * Constructs the link of 'simplex' with points coordinates.
-   */
-  Link_complex link(Edge_handle edge) const {
-    return Link_complex(*this, edge);
-  }
-
-  /**
-   * Constructs the link of 'simplex' with points coordinates.
-   */
-  Link_complex link(const Simplex_handle& simplex) const {
-    return Link_complex(*this, simplex);
-  }
-
-  /**
-   * @brief Compute the local vertices of 's' in the current complex
-   * If one of them is not present in the complex then the return value is uninitialized.
-   *
-   *
-   */
-  // xxx rename get_address et place un using dans sub_complex
-
-  boost::optional<Simplex_handle> get_simplex_address(
-                                                      const Root_simplex_handle& s) const {
-    boost::optional<Simplex_handle> res;
-
-    Simplex_handle s_address;
-    // Root_simplex_const_iterator i;
-    for (auto i = s.begin(); i != s.end(); ++i) {
-      boost::optional<Vertex_handle> address = get_address(*i);
-      if (!address)
-        return res;
-      else
-        s_address.add_vertex(*address);
-    }
-    res = s_address;
-    return res;
-  }
-
-  /**
-   * @brief returns a simplex with vertices which are the id of vertices of the
-   * argument.
-   */
-  Root_simplex_handle get_id(const Simplex_handle& local_simplex) const {
-    Root_simplex_handle global_simplex;
-    for (auto x = local_simplex.begin(); x != local_simplex.end(); ++x) {
-      global_simplex.add_vertex(get_id(*x));
-    }
-    return global_simplex;
-  }
-
-  /**
-   * @brief returns true iff the simplex s belongs to the simplicial
-   * complex.
-   */
-  virtual bool contains(const Simplex_handle & s) const {
-    if (s.dimension() == -1) {
-      return false;
-    } else if (s.dimension() == 0) {
-      return contains_vertex(s.first_vertex());
-    } else {
-      return (contains_edges(s) && !blocks(s));
-    }
-  }
-
-  /*
-   * @brief returnrs true iff the complex is empty.
-   */
-  bool empty() const {
-    return num_vertices() == 0;
-  }
-
-  /*
-   * @brief returns the number of vertices in the complex.
-   */
-  int num_vertices() const {
-    // remark boost::num_vertices(skeleton) counts deactivated vertices
-    return num_vertices_;
-  }
-
-  /*
-   * @brief returns the number of edges in the complex.
-   * @details currently in O(n)
-   */
-  // todo cache the value
-
-  int num_edges() const {
-    return boost::num_edges(skeleton);
-  }
-
-  /*
-   * @brief returns the number of blockers in the complex.
-   */
-  int num_blockers() const {
-    return num_blockers_;
-  }
-
-  /*
-   * @brief returns true iff the graph of the 1-skeleton of the complex is complete.
-   */
-  bool complete() const {
-    return (num_vertices() * (num_vertices() - 1)) / 2 == num_edges();
-  }
-
-  /**
-   * @brief returns the number of connected components in the graph of the 1-skeleton.
-   */
-  int num_connected_components() const {
-    int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices();
-    std::vector<int> component(skeleton.vertex_set().size());
-    return boost::connected_components(this->skeleton, &component[0])
-        - num_vert_collapsed;
-  }
-
-  /**
-   * @brief %Test if the complex is a cone.
-   * @details Runs in O(n) where n is the number of vertices.
-   */
-  bool is_cone() const {
-    if (num_vertices() == 0)
-      return false;
-    if (num_vertices() == 1)
-      return true;
-    for (auto vi : vertex_range()) {
-      // xxx todo faire une methode bool is_in_blocker(Vertex_handle)
-      if (blocker_map_.find(vi) == blocker_map_.end()) {
-        // no blocker passes through the vertex, we just need to
-        // check if the current vertex is linked to all others vertices of the complex
-        if (degree_[vi.vertex] == num_vertices() - 1)
-          return true;
-      }
-    }
-    return false;
-  }
-
-  //@}
+	template<class ComplexType> friend class Complex_vertex_iterator;
+	template<class ComplexType> friend class Complex_neighbors_vertices_iterator;
+	template<class ComplexType> friend class Complex_edge_iterator;
+	template<class ComplexType> friend class Complex_edge_around_vertex_iterator;
+
+	template<class ComplexType> friend class Skeleton_blocker_link_complex;
+	template<class ComplexType> friend class Skeleton_blocker_link_superior;
+	template<class ComplexType> friend class Skeleton_blocker_sub_complex;
+
+public:
+	/**
+	 * @brief The type of stored vertex node, specified by the template SkeletonBlockerDS
+	 */
+	typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex;
+
+	/**
+	 * @brief The type of stored  edge node, specified by the template SkeletonBlockerDS
+	 */
+	typedef typename SkeletonBlockerDS::Graph_edge Graph_edge;
+
+	typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle;
+
+	/**
+	 * @brief The type of an handle to a vertex of the complex.
+	 */
+	typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle;
+	typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle;
+
+	/**
+	 * @brief A ordered set of integers that represents a simplex.
+	 */
+	typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle;
+	typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle;
+
+	/**
+	 * @brief Handle to a blocker of the complex.
+	 */
+	typedef Simplex_handle* Blocker_handle;
+
+	typedef typename Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_iterator;
+	typedef typename Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator;
+
+protected:
+	typedef typename boost::adjacency_list<boost::setS,  // edges
+			boost::vecS,  // vertices
+			boost::undirectedS, Graph_vertex, Graph_edge> Graph;
+	// todo/remark : edges are not sorted, it heavily penalizes computation for SuperiorLink
+	// (eg Link with greater vertices)
+	// that burdens simplex iteration / complex initialization via list of simplices.
+	// to avoid that, one should modify the graph by storing two lists of adjacency for every
+	// vertex, the one with superior and the one with lower vertices, that way there is
+	// no more extra cost for computation of SuperiorLink
+	typedef typename boost::graph_traits<Graph>::vertex_iterator boost_vertex_iterator;
+	typedef typename boost::graph_traits<Graph>::edge_iterator boost_edge_iterator;
+
+protected:
+	typedef typename boost::graph_traits<Graph>::adjacency_iterator boost_adjacency_iterator;
+
+public:
+	/**
+	 * @brief Handle to an edge of the complex.
+	 */
+	typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle;
+
+protected:
+	typedef std::multimap<Vertex_handle, Simplex_handle *> BlockerMap;
+	typedef typename std::multimap<Vertex_handle, Simplex_handle *>::value_type BlockerPair;
+	typedef typename std::multimap<Vertex_handle, Simplex_handle *>::iterator BlockerMapIterator;
+	typedef typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator BlockerMapConstIterator;
+
+protected:
+	int num_vertices_;
+	int num_blockers_;
+
+	typedef Skeleton_blocker_complex_visitor<Vertex_handle> Visitor;
+	// typedef Visitor* Visitor_ptr;
+	Visitor* visitor;
+
+	/**
+	 * @details If 'x' is a Vertex_handle of a vertex in the complex then degree[x] = d is its degree.
+	 *
+	 * This quantity is updated when adding/removing edge.
+	 *
+	 * This is useful because the operation
+	 * list.size() is done in linear time.
+	 */
+	std::vector<boost_vertex_handle> degree_;
+	Graph skeleton; /** 1-skeleton of the simplicial complex. */
+
+	/** Each vertex can access to the blockers passing through it. */
+	BlockerMap blocker_map_;
+
+public:
+	/////////////////////////////////////////////////////////////////////////////
+	/** @name Constructors, Destructors
+	 */
+	//@{
+
+	/**
+	 *@brief constructs a simplicial complex with a given number of vertices and a visitor.
+	 */
+	explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL)
+	: visitor(visitor_) {
+		clear();
+		for (int i = 0; i < num_vertices_; ++i) {
+			add_vertex();
+		}
+	}
+
+private:
+	typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie;
+
+	template<typename SimpleHandleOutputIterator>
+	/**
+	 * return a vector of simplex trie for every vertex
+	 */
+	std::vector<STrie*> compute_cofaces(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) {
+		std::vector<STrie*> cofaces(num_vertices(), 0);
+		for (auto i = 0u; i < num_vertices(); ++i)
+			cofaces[i] = new STrie(Vertex_handle(i));
+		for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) {
+			if (s_it->dimension() >= 1)
+				cofaces[s_it->first_vertex()]->add_simplex(*s_it);
+		}
+		return cofaces;
+	}
+
+
+public:
+	/**
+	 * @brief Constructor with a list of simplices.
+	 * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not).
+	 */
+	template<typename SimpleHandleOutputIterator>
+	Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, bool is_flag_complex = false, Visitor* visitor_ = NULL)
+	: num_vertices_(0),
+	  num_blockers_(0),
+	  visitor(visitor_) {
+		std::vector<std::pair<Vertex_handle, Vertex_handle>> edges;
+		// first pass to add vertices and edges
+		int num_vertex = -1;
+		for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) {
+			if (s_it->dimension() == 0) num_vertex = (std::max)(num_vertex, s_it->first_vertex().vertex);
+			if (s_it->dimension() == 1) edges.emplace_back(s_it->first_vertex(), s_it->last_vertex());
+		}
+		while (num_vertex-- >= 0) add_vertex();
+
+		for (const auto& e : edges)
+			add_edge(e.first, e.second);
+
+		if (!is_flag_complex) {
+			// need to compute blockers
+			std::vector<STrie*> cofaces(compute_cofaces(simplex_begin, simplex_end));
+			std::vector<Simplex_handle> max_faces;
+
+			for (auto i = 0u; i < num_vertices(); ++i) {
+				auto max_faces_i = cofaces[i]->maximal_faces();
+				max_faces.insert(max_faces.end(), max_faces_i.begin(), max_faces_i.end());
+			}
+
+			// all time here
+
+			// for every maximal face s, it picks its first vertex v and check for all nv \in N(v)
+			// if s union nv is in the complex, if not it is a blocker.
+			for (auto max_face : max_faces) {
+				if (max_face.dimension() > 0) {
+					auto first_v = max_face.first_vertex();
+					for (auto nv : vertex_range(first_v)) {
+						if (!max_face.contains(nv) && max_face.first_vertex() < nv) {
+							// check that all edges in vertices(max_face)\cup nv are here
+							// since max_face is a simplex, we only need to check that edges
+							// (x nv) where x \in max_face are present
+							bool all_edges_here = true;
+							for (auto x : max_face)
+								if (!contains_edge(x, nv)) {
+									all_edges_here = false;
+									break;
+								}
+							if (all_edges_here) {  // eg  this->contains(max_face)
+								max_face.add_vertex(nv);
+								if (!cofaces[first_v]->contains(max_face)) {
+									// if there exists a blocker included in max_face, we remove it
+									// as it is not a minimum missing face
+									// the other alternative would be to check to all properfaces
+									// are in the complex before adding a blocker but that
+									// would be more expensive if there are few blockers
+									std::vector<Blocker_handle> blockers_to_remove;
+									bool is_blocker = true;
+									for (auto b : blocker_range(first_v))
+										if (b->contains(max_face)){
+											blockers_to_remove.push_back(b);
+										}
+										else
+											if (max_face.contains(*b)) {
+												is_blocker = false;
+												break;
+											}
+									for (auto b : blockers_to_remove)
+										this->delete_blocker(b);
+
+									if(is_blocker)
+										add_blocker(max_face);
+								}
+								max_face.remove_vertex(nv);
+							}
+						}
+					}
+				}
+			}
+			for (auto i = 0u; i < num_vertices(); ++i)
+				delete(cofaces[i]);
+		}
+	}
+
+	// We cannot use the default copy constructor since we need
+	// to make a copy of each of the blockers
+
+	Skeleton_blocker_complex(const Skeleton_blocker_complex& copy) {
+		visitor = NULL;
+		degree_ = copy.degree_;
+		skeleton = Graph(copy.skeleton);
+		num_vertices_ = copy.num_vertices_;
+
+		num_blockers_ = 0;
+		// we copy the blockers
+		for (auto blocker : copy.const_blocker_range()) {
+			add_blocker(*blocker);
+		}
+	}
+
+	/**
+	 */
+	Skeleton_blocker_complex& operator=(const Skeleton_blocker_complex& copy) {
+		clear();
+		visitor = NULL;
+		degree_ = copy.degree_;
+		skeleton = Graph(copy.skeleton);
+		num_vertices_ = copy.num_vertices_;
+
+		num_blockers_ = 0;
+		// we copy the blockers
+		for (auto blocker : copy.const_blocker_range())
+			add_blocker(*blocker);
+		return *this;
+	}
+
+
+	/**
+	 * return true if both complexes have the same simplices.
+	 */
+	bool operator==(const Skeleton_blocker_complex& other) const{
+		if(other.num_vertices() != num_vertices()) return false;
+		if(other.num_edges() != num_edges()) return false;
+		if(other.num_blockers() != num_blockers()) return false;
+
+		for(auto v : vertex_range())
+			if(!other.contains_vertex(v)) return false;
+
+		for(auto e : edge_range())
+			if(!other.contains_edge(first_vertex(e),second_vertex(e))) return false;
+
+		for(const auto b : const_blocker_range())
+			if(!other.contains_blocker(*b)) return false;
+
+		return true;
+	}
+
+	bool operator!=(const Skeleton_blocker_complex& other) const{
+		return !(*this == other);
+	}
+
+	/**
+	 * The destructor delete all blockers allocated.
+	 */
+	virtual ~Skeleton_blocker_complex() {
+		clear();
+	}
+
+	/**
+	 * @details Clears the simplicial complex. After a call to this function,
+	 * blockers are destroyed. The 1-skeleton and the set of blockers
+	 * are both empty.
+	 */
+	virtual void clear() {
+		// xxx for now the responsabilty of freeing the visitor is for
+		// the user
+		visitor = NULL;
+
+		degree_.clear();
+		num_vertices_ = 0;
+
+		remove_blockers();
+
+		skeleton.clear();
+	}
+
+	/**
+	 *@brief allows to change the visitor.
+	 */
+	void set_visitor(Visitor* other_visitor) {
+		visitor = other_visitor;
+	}
+
+	//@}
+
+	/////////////////////////////////////////////////////////////////////////////
+	/** @name Vertices operations
+	 */
+	//@{
+public:
+	/**
+	 * @brief Return a local Vertex_handle of a vertex given a global one.
+	 * @remark Assume that the vertex is present in the complex.
+	 */
+	Vertex_handle operator[](Root_vertex_handle global) const {
+		auto local(get_address(global));
+		assert(local);
+		return *local;
+	}
+
+	/**
+	 * @brief Return the vertex node associated to local Vertex_handle.
+	 * @remark Assume that the vertex is present in the complex.
+	 */
+	Graph_vertex& operator[](Vertex_handle address) {
+		assert(
+				0 <= address.vertex && address.vertex < boost::num_vertices(skeleton));
+		return skeleton[address.vertex];
+	}
+
+	/**
+	 * @brief Return the vertex node associated to local Vertex_handle.
+	 * @remark Assume that the vertex is present in the complex.
+	 */
+	const Graph_vertex& operator[](Vertex_handle address) const {
+		assert(
+				0 <= address.vertex && address.vertex < boost::num_vertices(skeleton));
+		return skeleton[address.vertex];
+	}
+
+	/**
+	 * @brief Adds a vertex to the simplicial complex and returns its Vertex_handle.
+	 */
+	Vertex_handle add_vertex() {
+		Vertex_handle address(boost::add_vertex(skeleton));
+		num_vertices_++;
+		(*this)[address].activate();
+		// safe since we now that we are in the root complex and the field 'address' and 'id'
+		// are identical for every vertices
+		(*this)[address].set_id(Root_vertex_handle(address.vertex));
+		degree_.push_back(0);
+		if (visitor)
+			visitor->on_add_vertex(address);
+		return address;
+	}
+
+	/**
+	 * @brief Remove a vertex from the simplicial complex
+	 * @remark It just deactivates the vertex with a boolean flag but does not
+	 * remove it from vertices from complexity issues.
+	 */
+	void remove_vertex(Vertex_handle address) {
+		assert(contains_vertex(address));
+		// We remove b
+		boost::clear_vertex(address.vertex, skeleton);
+		(*this)[address].deactivate();
+		num_vertices_--;
+		degree_[address.vertex] = -1;
+		if (visitor)
+			visitor->on_remove_vertex(address);
+	}
+
+	/**
+	 */
+	bool contains_vertex(Vertex_handle u) const {
+		if (u.vertex < 0 || u.vertex >= boost::num_vertices(skeleton))
+			return false;
+		return (*this)[u].is_active();
+	}
+
+	/**
+	 */
+	bool contains_vertex(Root_vertex_handle u) const {
+		boost::optional<Vertex_handle> address = get_address(u);
+		return address && (*this)[*address].is_active();
+	}
+
+	/**
+	 * @return true iff the simplicial complex contains all vertices
+	 * of simplex sigma
+	 */
+	bool contains_vertices(const Simplex_handle & sigma) const {
+		for (auto vertex : sigma)
+			if (!contains_vertex(vertex))
+				return false;
+		return true;
+	}
+
+	/**
+	 * @brief Given an Id return the address of the vertex having this Id in the complex.
+	 * @remark For a simplicial complex, the address is the id but it may not be the case for a SubComplex.
+	 */
+	virtual boost::optional<Vertex_handle> get_address(
+			Root_vertex_handle id) const {
+		boost::optional<Vertex_handle> res;
+		if (id.vertex < boost::num_vertices(skeleton))
+			res = Vertex_handle(id.vertex);  // xxx
+		return res;
+	}
+
+	/**
+	 * return the id of a vertex of adress local present in the graph
+	 */
+	Root_vertex_handle get_id(Vertex_handle local) const {
+		assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton));
+		return (*this)[local].get_id();
+	}
+
+	/**
+	 * @brief Convert an address of a vertex of a complex to the address in
+	 * the current complex.
+	 * @details
+	 * If the current complex is a sub (or sup) complex of 'other', it converts
+	 * the address of a vertex v expressed in 'other' to the address of the vertex
+	 * v in the current one.
+	 * @remark this methods uses Root_vertex_handle to identify the vertex and
+	 * assumes the vertex is present in the current complex.
+	 */
+	Vertex_handle convert_handle_from_another_complex(
+			const Skeleton_blocker_complex& other, Vertex_handle vh_in_other) const {
+		auto vh_in_current_complex = get_address(other.get_id(vh_in_other));
+		assert(vh_in_current_complex);
+		return *vh_in_current_complex;
+	}
+
+	/**
+	 * @brief return the graph degree of a vertex.
+	 */
+	int degree(Vertex_handle local) const {
+		assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton));
+		return degree_[local.vertex];
+	}
+
+	//@}
+
+	/////////////////////////////////////////////////////////////////////////////
+	/** @name Edges operations
+	 */
+	//@{
+public:
+	/**
+	 * @brief return an edge handle if the two vertices forms
+	 * an edge in the complex
+	 */
+	boost::optional<Edge_handle> operator[](
+			const std::pair<Vertex_handle, Vertex_handle>& ab) const {
+		boost::optional<Edge_handle> res;
+		std::pair<Edge_handle, bool> edge_pair(
+				boost::edge(ab.first.vertex, ab.second.vertex, skeleton));
+		if (edge_pair.second)
+			res = edge_pair.first;
+		return res;
+	}
+
+	/**
+	 * @brief returns the stored node associated to an edge
+	 */
+	Graph_edge& operator[](Edge_handle edge_handle) {
+		return skeleton[edge_handle];
+	}
+
+	/**
+	 * @brief returns the stored node associated to an edge
+	 */
+	const Graph_edge& operator[](Edge_handle edge_handle) const {
+		return skeleton[edge_handle];
+	}
+
+	/**
+	 * @brief returns the first vertex of an edge
+	 * @details it assumes that the edge is present in the complex
+	 */
+	Vertex_handle first_vertex(Edge_handle edge_handle) const {
+		return static_cast<Vertex_handle> (source(edge_handle, skeleton));
+	}
+
+	/**
+	 * @brief returns the first vertex of an edge
+	 * @details it assumes that the edge is present in the complex
+	 */
+	Vertex_handle second_vertex(Edge_handle edge_handle) const {
+		return static_cast<Vertex_handle> (target(edge_handle, skeleton));
+	}
+
+	/**
+	 * @brief returns the simplex made with the two vertices of the edge
+	 * @details it assumes that the edge is present in the complex
+
+	 */
+	Simplex_handle get_vertices(Edge_handle edge_handle) const {
+		auto edge((*this)[edge_handle]);
+		return Simplex_handle((*this)[edge.first()], (*this)[edge.second()]);
+	}
+
+	/**
+	 * @brief Adds an edge between vertices a and b and all its cofaces.
+	 */
+	Edge_handle add_edge(Vertex_handle a, Vertex_handle b) {
+		assert(contains_vertex(a) && contains_vertex(b));
+		assert(a != b);
+
+		auto edge_handle((*this)[std::make_pair(a, b)]);
+		// std::pair<Edge_handle,bool> pair_descr_bool = (*this)[std::make_pair(a,b)];
+		// Edge_handle edge_descr;
+		// bool edge_present = pair_descr_bool.second;
+		if (!edge_handle) {
+			edge_handle = boost::add_edge(a.vertex, b.vertex, skeleton).first;
+			(*this)[*edge_handle].setId(get_id(a), get_id(b));
+			degree_[a.vertex]++;
+			degree_[b.vertex]++;
+			if (visitor)
+				visitor->on_add_edge(a, b);
+		}
+		return *edge_handle;
+	}
+
+	/**
+	 * @brief Adds all edges and their cofaces of a simplex to the simplicial complex.
+	 */
+	void add_edges(const Simplex_handle & sigma) {
+		Simplex_handle_iterator i, j;
+		for (i = sigma.begin(); i != sigma.end(); ++i)
+			for (j = i, j++; j != sigma.end(); ++j)
+				add_edge(*i, *j);
+	}
+
+	/**
+	 * @brief Removes an edge from the simplicial complex and all its cofaces.
+	 * @details returns the former Edge_handle representing the edge
+	 */
+	virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b) {
+		bool found;
+		Edge_handle edge;
+		tie(edge, found) = boost::edge(a.vertex, b.vertex, skeleton);
+		if (found) {
+			if (visitor)
+				visitor->on_remove_edge(a, b);
+			// if (heapCollapse.Contains(edge)) heapCollapse.Delete(edge);
+			boost::remove_edge(a.vertex, b.vertex, skeleton);
+			degree_[a.vertex]--;
+			degree_[b.vertex]--;
+		}
+		return edge;
+	}
+
+	/**
+	 * @brief Removes edge and its cofaces from the simplicial complex.
+	 */
+	void remove_edge(Edge_handle edge) {
+		assert(contains_vertex(first_vertex(edge)));
+		assert(contains_vertex(second_vertex(edge)));
+		remove_edge(first_vertex(edge), second_vertex(edge));
+	}
+
+	/**
+	 * @brief The complex is reduced to its set of vertices.
+	 * All the edges and blockers are removed.
+	 */
+	void keep_only_vertices() {
+		remove_blockers();
+
+		for (auto u : vertex_range()) {
+			while (this->degree(u) > 0) {
+				Vertex_handle v(*(adjacent_vertices(u.vertex, this->skeleton).first));
+				this->remove_edge(u, v);
+			}
+		}
+	}
+
+	/**
+	 * @return true iff the simplicial complex contains an edge between
+	 * vertices a and b
+	 */
+	bool contains_edge(Vertex_handle a, Vertex_handle b) const {
+		// if (a.vertex<0 || b.vertex <0) return false;
+		return boost::edge(a.vertex, b.vertex, skeleton).second;
+	}
+
+	/**
+	 * @return true iff the simplicial complex contains all vertices
+	 * and all edges of simplex sigma
+	 */
+	bool contains_edges(const Simplex_handle & sigma) const {
+		for (auto i = sigma.begin(); i != sigma.end(); ++i) {
+			if (!contains_vertex(*i))
+				return false;
+			for (auto j = i; ++j != sigma.end();) {
+				if (!contains_edge(*i, *j))
+					return false;
+			}
+		}
+		return true;
+	}
+	//@}
+
+	/////////////////////////////////////////////////////////////////////////////
+	/** @name Blockers operations
+	 */
+	//@{
+
+	/**
+	 * @brief Adds the simplex to the set of blockers and
+	 * returns a Blocker_handle toward it if was not present before and 0 otherwise.
+	 */
+	Blocker_handle add_blocker(const Simplex_handle& blocker) {
+		assert(blocker.dimension() > 1);
+		if (contains_blocker(blocker)) {
+			// std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl;
+			return 0;
+		} else {
+			if (visitor)
+				visitor->on_add_blocker(blocker);
+			Blocker_handle blocker_pt = new Simplex_handle(blocker);
+			num_blockers_++;
+			auto vertex = blocker_pt->begin();
+			while (vertex != blocker_pt->end()) {
+				blocker_map_.insert(BlockerPair(*vertex, blocker_pt));
+				++vertex;
+			}
+			return blocker_pt;
+		}
+	}
+
+protected:
+	/**
+	 * @brief Adds the simplex to the set of blockers
+	 */
+	void add_blocker(Blocker_handle blocker) {
+		if (contains_blocker(*blocker)) {
+			// std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl;
+			return;
+		} else {
+			if (visitor)
+				visitor->on_add_blocker(*blocker);
+			num_blockers_++;
+			auto vertex = blocker->begin();
+			while (vertex != blocker->end()) {
+				blocker_map_.insert(BlockerPair(*vertex, blocker));
+				++vertex;
+			}
+		}
+	}
+
+protected:
+	/**
+	 * Removes sigma from the blocker map of vertex v
+	 */
+	void remove_blocker(const Blocker_handle sigma, Vertex_handle v) {
+		Complex_blocker_around_vertex_iterator blocker;
+		for (blocker = blocker_range(v).begin(); blocker != blocker_range(v).end();
+				++blocker) {
+			if (*blocker == sigma)
+				break;
+		}
+		if (*blocker != sigma) {
+			std::cerr
+			<< "bug ((*blocker).second == sigma) ie try to remove a blocker not present\n";
+			assert(false);
+		} else {
+			blocker_map_.erase(blocker.current_position());
+		}
+	}
+
+public:
+	/**
+	 * @brief Removes the simplex from the set of blockers.
+	 * @remark sigma has to belongs to the set of blockers
+	 */
+	void remove_blocker(const Blocker_handle sigma) {
+		for (auto vertex : *sigma)
+			remove_blocker(sigma, vertex);
+		num_blockers_--;
+	}
+
+	/**
+	 * @brief Remove all blockers, in other words, it expand the simplicial
+	 * complex to the smallest flag complex that contains it.
+	 */
+	void remove_blockers() {
+		// Desallocate the blockers
+		while (!blocker_map_.empty()) {
+			delete_blocker(blocker_map_.begin()->second);
+		}
+		num_blockers_ = 0;
+		blocker_map_.clear();
+	}
+
+protected:
+	/**
+	 * Removes the simplex sigma from the set of blockers.
+	 * sigma has to belongs to the set of blockers
+	 *
+	 * @remark contrarily to delete_blockers does not call the destructor
+	 */
+	void remove_blocker(const Simplex_handle& sigma) {
+		assert(contains_blocker(sigma));
+		for (auto vertex : sigma)
+			remove_blocker(sigma, vertex);
+		num_blockers_--;
+	}
+
+public:
+	/**
+	 * Removes the simplex s from the set of blockers
+	 * and desallocate s.
+	 */
+	void delete_blocker(Blocker_handle sigma) {
+		if (visitor)
+			visitor->on_delete_blocker(sigma);
+		remove_blocker(sigma);
+		delete sigma;
+	}
+
+	/**
+	 * @return true iff s is a blocker of the simplicial complex
+	 */
+	bool contains_blocker(const Blocker_handle s) const {
+		if (s->dimension() < 2)
+			return false;
+
+		Vertex_handle a = s->first_vertex();
+
+		for (const auto blocker : const_blocker_range(a)) {
+			if (s == *blocker)
+				return true;
+		}
+		return false;
+	}
+
+	/**
+	 * @return true iff s is a blocker of the simplicial complex
+	 */
+	bool contains_blocker(const Simplex_handle & s) const {
+		if (s.dimension() < 2)
+			return false;
+
+		Vertex_handle a = s.first_vertex();
+
+		for (auto blocker : const_blocker_range(a)) {
+			if (s == *blocker)
+				return true;
+		}
+		return false;
+	}
+
+private:
+	/**
+	 * @return true iff a blocker of the simplicial complex
+	 * is a face of sigma.
+	 */
+	bool blocks(const Simplex_handle & sigma) const {
+		for (auto blocker : const_blocker_range()) {
+			if (sigma.contains(*blocker))
+				return true;
+		}
+		return false;
+	}
+
+	//@}
+
+protected:
+	/**
+	 * @details Adds to simplex the neighbours of v e.g. \f$ n \leftarrow n \cup N(v) \f$.
+	 * If keep_only_superior is true then only vertices that are greater than v are added.
+	 */
+	virtual void add_neighbours(Vertex_handle v, Simplex_handle & n,
+			bool keep_only_superior = false) const {
+		boost_adjacency_iterator ai, ai_end;
+		for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end;
+				++ai) {
+			if (keep_only_superior) {
+				if (*ai > v.vertex) {
+					n.add_vertex(Vertex_handle(*ai));
+				}
+			} else {
+				n.add_vertex(Vertex_handle(*ai));
+			}
+		}
+	}
+
+	/**
+	 * @details Add to simplex res all vertices which are
+	 * neighbours of alpha: ie \f$ res \leftarrow res \cup N(alpha) \f$.
+	 *
+	 * If 'keep_only_superior' is true then only vertices that are greater than alpha are added.
+	 * todo revoir
+	 *
+	 */
+	virtual void add_neighbours(const Simplex_handle &alpha, Simplex_handle & res,
+			bool keep_only_superior = false) const {
+		res.clear();
+		auto alpha_vertex = alpha.begin();
+		add_neighbours(*alpha_vertex, res, keep_only_superior);
+		for (alpha_vertex = (alpha.begin())++; alpha_vertex != alpha.end();
+				++alpha_vertex)
+			keep_neighbours(*alpha_vertex, res, keep_only_superior);
+	}
+
+	/**
+	 * @details Remove from simplex n all vertices which are
+	 * not neighbours of v e.g. \f$ res \leftarrow res \cap N(v) \f$.
+	 * If 'keep_only_superior' is true then only vertices that are greater than v are keeped.
+	 */
+	virtual void keep_neighbours(Vertex_handle v, Simplex_handle& res,
+			bool keep_only_superior = false) const {
+		Simplex_handle nv;
+		add_neighbours(v, nv, keep_only_superior);
+		res.intersection(nv);
+	}
+
+	/**
+	 * @details Remove from simplex all vertices which are
+	 * neighbours of v eg \f$ res \leftarrow res \setminus N(v) \f$.
+	 * If 'keep_only_superior' is true then only vertices that are greater than v are added.
+	 */
+	virtual void remove_neighbours(Vertex_handle v, Simplex_handle & res,
+			bool keep_only_superior = false) const {
+		Simplex_handle nv;
+		add_neighbours(v, nv, keep_only_superior);
+		res.difference(nv);
+	}
+
+public:
+	typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex> Link_complex;
+
+	/**
+	 * Constructs the link of 'simplex' with points coordinates.
+	 */
+	Link_complex link(Vertex_handle v) const {
+		return Link_complex(*this, Simplex_handle(v));
+	}
+
+	/**
+	 * Constructs the link of 'simplex' with points coordinates.
+	 */
+	Link_complex link(Edge_handle edge) const {
+		return Link_complex(*this, edge);
+	}
+
+	/**
+	 * Constructs the link of 'simplex' with points coordinates.
+	 */
+	Link_complex link(const Simplex_handle& simplex) const {
+		return Link_complex(*this, simplex);
+	}
+
+	/**
+	 * @brief Compute the local vertices of 's' in the current complex
+	 * If one of them is not present in the complex then the return value is uninitialized.
+	 *
+	 *
+	 */
+	// xxx rename get_address et place un using dans sub_complex
+
+	boost::optional<Simplex_handle> get_simplex_address(
+			const Root_simplex_handle& s) const {
+		boost::optional<Simplex_handle> res;
+
+		Simplex_handle s_address;
+		// Root_simplex_const_iterator i;
+		for (auto i = s.begin(); i != s.end(); ++i) {
+			boost::optional<Vertex_handle> address = get_address(*i);
+			if (!address)
+				return res;
+			else
+				s_address.add_vertex(*address);
+		}
+		res = s_address;
+		return res;
+	}
+
+	/**
+	 * @brief returns a simplex with vertices which are the id of vertices of the
+	 * argument.
+	 */
+	Root_simplex_handle get_id(const Simplex_handle& local_simplex) const {
+		Root_simplex_handle global_simplex;
+		for (auto x = local_simplex.begin(); x != local_simplex.end(); ++x) {
+			global_simplex.add_vertex(get_id(*x));
+		}
+		return global_simplex;
+	}
+
+	/**
+	 * @brief returns true iff the simplex s belongs to the simplicial
+	 * complex.
+	 */
+	virtual bool contains(const Simplex_handle & s) const {
+		if (s.dimension() == -1) {
+			return false;
+		} else if (s.dimension() == 0) {
+			return contains_vertex(s.first_vertex());
+		} else {
+			return (contains_edges(s) && !blocks(s));
+		}
+	}
+
+	/*
+	 * @brief returnrs true iff the complex is empty.
+	 */
+	bool empty() const {
+		return num_vertices() == 0;
+	}
+
+	/*
+	 * @brief returns the number of vertices in the complex.
+	 */
+	int num_vertices() const {
+		// remark boost::num_vertices(skeleton) counts deactivated vertices
+		return num_vertices_;
+	}
+
+	/*
+	 * @brief returns the number of edges in the complex.
+	 * @details currently in O(n)
+	 */
+	// todo cache the value
+
+	int num_edges() const {
+		return boost::num_edges(skeleton);
+	}
+
+	int num_triangles() const{
+		auto triangles = triangle_range();
+		return std::distance(triangles.begin(),triangles.end());
+	}
+	/*
+	 * @brief returns the number of blockers in the complex.
+	 */
+	int num_blockers() const {
+		return num_blockers_;
+	}
+
+	/*
+	 * @brief returns true iff the graph of the 1-skeleton of the complex is complete.
+	 */
+	bool complete() const {
+		return (num_vertices() * (num_vertices() - 1)) / 2 == num_edges();
+	}
+
+	/**
+	 * @brief returns the number of connected components in the graph of the 1-skeleton.
+	 */
+	int num_connected_components() const {
+		int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices();
+		std::vector<int> component(skeleton.vertex_set().size());
+		return boost::connected_components(this->skeleton, &component[0])
+		- num_vert_collapsed;
+	}
+
+	/**
+	 * @brief %Test if the complex is a cone.
+	 * @details Runs in O(n) where n is the number of vertices.
+	 */
+	bool is_cone() const {
+		if (num_vertices() == 0)
+			return false;
+		if (num_vertices() == 1)
+			return true;
+		for (auto vi : vertex_range()) {
+			// xxx todo faire une methode bool is_in_blocker(Vertex_handle)
+			if (blocker_map_.find(vi) == blocker_map_.end()) {
+				// no blocker passes through the vertex, we just need to
+				// check if the current vertex is linked to all others vertices of the complex
+				if (degree_[vi.vertex] == num_vertices() - 1)
+					return true;
+			}
+		}
+		return false;
+	}
+
+	//@}
 	/** @Simplification operations
 	 */
 	//@{
@@ -1244,279 +1255,279 @@ private:
 	//@}
 
 public:
-  /////////////////////////////////////////////////////////////////////////////
-  /** @name Vertex iterators
-   */
-  //@{
-  typedef Complex_vertex_iterator<Skeleton_blocker_complex> CVI;  // todo rename
-
-  //
-  // @brief Range over the vertices of the simplicial complex.
-  // Methods .begin() and .end() return a Complex_vertex_iterator.
-  //
-  typedef boost::iterator_range<
-  Complex_vertex_iterator<Skeleton_blocker_complex> > Complex_vertex_range;
-
-  /**
-   * @brief Returns a Complex_vertex_range over all vertices of the complex
-   */
-  Complex_vertex_range vertex_range() const {
-    auto begin = Complex_vertex_iterator<Skeleton_blocker_complex>(this);
-    auto end = Complex_vertex_iterator<Skeleton_blocker_complex>(this, 0);
-    return Complex_vertex_range(begin, end);
-  }
-
-  typedef boost::iterator_range<
-  Complex_neighbors_vertices_iterator<Skeleton_blocker_complex> > Complex_neighbors_vertices_range;
-
-  /**
-   * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v
-   */
-  Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const {
-    auto begin = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(
-                                                                               this, v);
-    auto end = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(
-                                                                             this, v, 0);
-    return Complex_neighbors_vertices_range(begin, end);
-  }
-
-  //@}
-
-  /** @name Edge iterators
-   */
-  //@{
-
-  typedef boost::iterator_range<
-  Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>> Complex_edge_range;
-
-  /**
-   * @brief Returns a Complex_edge_range over all edges of the simplicial complex
-   */
-  Complex_edge_range edge_range() const {
-    auto begin = Complex_edge_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this);
-    auto end = Complex_edge_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, 0);
-    return Complex_edge_range(begin, end);
-  }
-
-  typedef boost::iterator_range <Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>>
-  Complex_edge_around_vertex_range;
-
-  /**
-   * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes
-   * through 'v'
-   */
-  Complex_edge_around_vertex_range edge_range(Vertex_handle v) const {
-    auto begin = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, v);
-    auto end = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, v, 0);
-    return Complex_edge_around_vertex_range(begin, end);
-  }
-
-  //@}
-
-  /** @name Triangles iterators
-   */
-  //@{
- private:
-  typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link;
-  typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link;
-
- public:
-  typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex, Superior_link> Superior_triangle_around_vertex_iterator;
-
-  typedef boost::iterator_range < Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link> > Complex_triangle_around_vertex_range;
-
-  /**
-   * @brief Range over triangles around a vertex of the simplicial complex.
-   * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
-   *
-   */
-  Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const {
-    auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v);
-    auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v, 0);
-    return Complex_triangle_around_vertex_range(begin, end);
-  }
-
-  typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range;
-
-  /**
-   * @brief Range over triangles of the simplicial complex.
-   * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
-   *
-   */
-  Complex_triangle_range triangle_range() const {
-    auto end = Triangle_iterator<Skeleton_blocker_complex>(this, 0);
-    if (empty()) {
-      return Complex_triangle_range(end, end);
-    } else {
-      auto begin = Triangle_iterator<Skeleton_blocker_complex>(this);
-      return Complex_triangle_range(begin, end);
-    }
-  }
-
-  //@}
-
-  /** @name Simplices iterators
-   */
-  //@{
-  typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex, Link> Complex_simplex_around_vertex_iterator;
-
-  /**
-   * @brief Range over the simplices of the simplicial complex around a vertex.
-   * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator.
-   */
-  typedef boost::iterator_range < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range;
-
-  /**
-   * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex
-   */
-  Complex_simplex_around_vertex_range simplex_range(Vertex_handle v) const {
-    assert(contains_vertex(v));
-    return Complex_simplex_around_vertex_range(
-                                               Complex_simplex_around_vertex_iterator(this, v),
-                                               Complex_simplex_around_vertex_iterator(this, v, true));
-  }
-
-  // typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator;
-  typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator;
-
-  typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range;
-
-  /**
-   * @brief Returns a Complex_simplex_range over all the simplices of the complex
-   */
-  Complex_simplex_range simplex_range() const {
-    Complex_simplex_iterator end(this, true);
-    if (empty()) {
-      return Complex_simplex_range(end, end);
-    } else {
-      Complex_simplex_iterator begin(this);
-      return Complex_simplex_range(begin, end);
-    }
-  }
-
-  //@}
-
-  /** @name Blockers iterators
-   */
-  //@{
- private:
-  /**
-   * @brief Iterator over the blockers adjacent to a vertex
-   */
-  typedef Blocker_iterator_around_vertex_internal<
-  typename std::multimap<Vertex_handle, Simplex_handle *>::iterator,
-  Blocker_handle>
-  Complex_blocker_around_vertex_iterator;
-
-  /**
-   * @brief Iterator over (constant) blockers adjacent to a vertex
-   */
-  typedef Blocker_iterator_around_vertex_internal<
-  typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator,
-  const Blocker_handle>
-  Const_complex_blocker_around_vertex_iterator;
-
-  typedef boost::iterator_range <Complex_blocker_around_vertex_iterator> Complex_blocker_around_vertex_range;
-  typedef boost::iterator_range <Const_complex_blocker_around_vertex_iterator> Const_complex_blocker_around_vertex_range;
-
- public:
-  /**
-   * @brief Returns a range of the blockers of the complex passing through a vertex
-   */
-  Complex_blocker_around_vertex_range blocker_range(Vertex_handle v) {
-    auto begin = Complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v));
-    auto end = Complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v));
-    return Complex_blocker_around_vertex_range(begin, end);
-  }
-
-  /**
-   * @brief Returns a range of the blockers of the complex passing through a vertex
-   */
-  Const_complex_blocker_around_vertex_range const_blocker_range(Vertex_handle v) const {
-    auto begin = Const_complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v));
-    auto end = Const_complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v));
-    return Const_complex_blocker_around_vertex_range(begin, end);
-  }
-
- private:
-  /**
-   * @brief Iterator over the blockers.
-   */
-  typedef Blocker_iterator_internal<
-  typename std::multimap<Vertex_handle, Simplex_handle *>::iterator,
-  Blocker_handle>
-  Complex_blocker_iterator;
-
-  /**
-   * @brief Iterator over the (constant) blockers.
-   */
-  typedef Blocker_iterator_internal<
-  typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator,
-  const Blocker_handle>
-  Const_complex_blocker_iterator;
-
-  typedef boost::iterator_range <Complex_blocker_iterator> Complex_blocker_range;
-  typedef boost::iterator_range <Const_complex_blocker_iterator> Const_complex_blocker_range;
-
- public:
-  /**
-   * @brief Returns a range of the blockers of the complex
-   */
-  Complex_blocker_range blocker_range() {
-    auto begin = Complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end());
-    auto end = Complex_blocker_iterator(blocker_map_.end(), blocker_map_.end());
-    return Complex_blocker_range(begin, end);
-  }
-
-  /**
-   * @brief Returns a range of the blockers of the complex
-   */
-  Const_complex_blocker_range const_blocker_range() const {
-    auto begin = Const_complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end());
-    auto end = Const_complex_blocker_iterator(blocker_map_.end(), blocker_map_.end());
-    return Const_complex_blocker_range(begin, end);
-  }
-
-  //@}
-
-  /////////////////////////////////////////////////////////////////////////////
-  /** @name Print and IO methods
-   */
-  //@{
- public:
-  std::string to_string() const {
-    std::ostringstream stream;
-    stream << num_vertices() << " vertices:\n" << vertices_to_string() << std::endl;
-    stream << num_edges() << " edges:\n" << edges_to_string() << std::endl;
-    stream << num_blockers() << " blockers:\n" << blockers_to_string() << std::endl;
-    return stream.str();
-  }
-
-  std::string vertices_to_string() const {
-    std::ostringstream stream;
-    for (auto vertex : vertex_range()) {
-      stream << "{" << (*this)[vertex].get_id() << "} ";
-    }
-    stream << std::endl;
-    return stream.str();
-  }
-
-  std::string edges_to_string() const {
-    std::ostringstream stream;
-    for (auto edge : edge_range())
-      stream << "{" << (*this)[edge].first() << "," << (*this)[edge].second() << "} ";
-    stream << std::endl;
-    return stream.str();
-  }
-
-  std::string blockers_to_string() const {
-    std::ostringstream stream;
-
-    for (auto b : const_blocker_range())
-      stream << *b << std::endl;
-    return stream.str();
-  }
-  //@}
+	/////////////////////////////////////////////////////////////////////////////
+	/** @name Vertex iterators
+	 */
+	//@{
+	typedef Complex_vertex_iterator<Skeleton_blocker_complex> CVI;  // todo rename
+
+	//
+	// @brief Range over the vertices of the simplicial complex.
+	// Methods .begin() and .end() return a Complex_vertex_iterator.
+	//
+	typedef boost::iterator_range<
+			Complex_vertex_iterator<Skeleton_blocker_complex> > Complex_vertex_range;
+
+	/**
+	 * @brief Returns a Complex_vertex_range over all vertices of the complex
+	 */
+	Complex_vertex_range vertex_range() const {
+		auto begin = Complex_vertex_iterator<Skeleton_blocker_complex>(this);
+		auto end = Complex_vertex_iterator<Skeleton_blocker_complex>(this, 0);
+		return Complex_vertex_range(begin, end);
+	}
+
+	typedef boost::iterator_range<
+			Complex_neighbors_vertices_iterator<Skeleton_blocker_complex> > Complex_neighbors_vertices_range;
+
+	/**
+	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v
+	 */
+	Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const {
+		auto begin = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(
+				this, v);
+		auto end = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(
+				this, v, 0);
+		return Complex_neighbors_vertices_range(begin, end);
+	}
+
+	//@}
+
+	/** @name Edge iterators
+	 */
+	//@{
+
+	typedef boost::iterator_range<
+			Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>> Complex_edge_range;
+
+	/**
+	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex
+	 */
+	Complex_edge_range edge_range() const {
+		auto begin = Complex_edge_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this);
+		auto end = Complex_edge_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, 0);
+		return Complex_edge_range(begin, end);
+	}
+
+	typedef boost::iterator_range <Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>>
+			Complex_edge_around_vertex_range;
+
+	/**
+	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes
+	 * through 'v'
+	 */
+	Complex_edge_around_vertex_range edge_range(Vertex_handle v) const {
+		auto begin = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, v);
+		auto end = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex < SkeletonBlockerDS >> (this, v, 0);
+		return Complex_edge_around_vertex_range(begin, end);
+	}
+
+	//@}
+
+	/** @name Triangles iterators
+	 */
+	//@{
+private:
+	typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link;
+	typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link;
+
+public:
+	typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex, Superior_link> Superior_triangle_around_vertex_iterator;
+
+	typedef boost::iterator_range < Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link> > Complex_triangle_around_vertex_range;
+
+	/**
+	 * @brief Range over triangles around a vertex of the simplicial complex.
+	 * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
+	 *
+	 */
+	Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const {
+		auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v);
+		auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v, 0);
+		return Complex_triangle_around_vertex_range(begin, end);
+	}
+
+	typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range;
+
+	/**
+	 * @brief Range over triangles of the simplicial complex.
+	 * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
+	 *
+	 */
+	Complex_triangle_range triangle_range() const {
+		auto end = Triangle_iterator<Skeleton_blocker_complex>(this, 0);
+		if (empty()) {
+			return Complex_triangle_range(end, end);
+		} else {
+			auto begin = Triangle_iterator<Skeleton_blocker_complex>(this);
+			return Complex_triangle_range(begin, end);
+		}
+	}
+
+	//@}
+
+	/** @name Simplices iterators
+	 */
+	//@{
+	typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex, Link> Complex_simplex_around_vertex_iterator;
+
+	/**
+	 * @brief Range over the simplices of the simplicial complex around a vertex.
+	 * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator.
+	 */
+	typedef boost::iterator_range < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range;
+
+	/**
+	 * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex
+	 */
+	Complex_simplex_around_vertex_range simplex_range(Vertex_handle v) const {
+		assert(contains_vertex(v));
+		return Complex_simplex_around_vertex_range(
+				Complex_simplex_around_vertex_iterator(this, v),
+				Complex_simplex_around_vertex_iterator(this, v, true));
+	}
+
+	// typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator;
+	typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator;
+
+	typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range;
+
+	/**
+	 * @brief Returns a Complex_simplex_range over all the simplices of the complex
+	 */
+	Complex_simplex_range simplex_range() const {
+		Complex_simplex_iterator end(this, true);
+		if (empty()) {
+			return Complex_simplex_range(end, end);
+		} else {
+			Complex_simplex_iterator begin(this);
+			return Complex_simplex_range(begin, end);
+		}
+	}
+
+	//@}
+
+	/** @name Blockers iterators
+	 */
+	//@{
+private:
+	/**
+	 * @brief Iterator over the blockers adjacent to a vertex
+	 */
+	typedef Blocker_iterator_around_vertex_internal<
+			typename std::multimap<Vertex_handle, Simplex_handle *>::iterator,
+			Blocker_handle>
+	Complex_blocker_around_vertex_iterator;
+
+	/**
+	 * @brief Iterator over (constant) blockers adjacent to a vertex
+	 */
+	typedef Blocker_iterator_around_vertex_internal<
+			typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator,
+			const Blocker_handle>
+	Const_complex_blocker_around_vertex_iterator;
+
+	typedef boost::iterator_range <Complex_blocker_around_vertex_iterator> Complex_blocker_around_vertex_range;
+	typedef boost::iterator_range <Const_complex_blocker_around_vertex_iterator> Const_complex_blocker_around_vertex_range;
+
+public:
+	/**
+	 * @brief Returns a range of the blockers of the complex passing through a vertex
+	 */
+	Complex_blocker_around_vertex_range blocker_range(Vertex_handle v) {
+		auto begin = Complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v));
+		auto end = Complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v));
+		return Complex_blocker_around_vertex_range(begin, end);
+	}
+
+	/**
+	 * @brief Returns a range of the blockers of the complex passing through a vertex
+	 */
+	Const_complex_blocker_around_vertex_range const_blocker_range(Vertex_handle v) const {
+		auto begin = Const_complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v));
+		auto end = Const_complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v));
+		return Const_complex_blocker_around_vertex_range(begin, end);
+	}
+
+private:
+	/**
+	 * @brief Iterator over the blockers.
+	 */
+	typedef Blocker_iterator_internal<
+			typename std::multimap<Vertex_handle, Simplex_handle *>::iterator,
+			Blocker_handle>
+	Complex_blocker_iterator;
+
+	/**
+	 * @brief Iterator over the (constant) blockers.
+	 */
+	typedef Blocker_iterator_internal<
+			typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator,
+			const Blocker_handle>
+	Const_complex_blocker_iterator;
+
+	typedef boost::iterator_range <Complex_blocker_iterator> Complex_blocker_range;
+	typedef boost::iterator_range <Const_complex_blocker_iterator> Const_complex_blocker_range;
+
+public:
+	/**
+	 * @brief Returns a range of the blockers of the complex
+	 */
+	Complex_blocker_range blocker_range() {
+		auto begin = Complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end());
+		auto end = Complex_blocker_iterator(blocker_map_.end(), blocker_map_.end());
+		return Complex_blocker_range(begin, end);
+	}
+
+	/**
+	 * @brief Returns a range of the blockers of the complex
+	 */
+	Const_complex_blocker_range const_blocker_range() const {
+		auto begin = Const_complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end());
+		auto end = Const_complex_blocker_iterator(blocker_map_.end(), blocker_map_.end());
+		return Const_complex_blocker_range(begin, end);
+	}
+
+	//@}
+
+	/////////////////////////////////////////////////////////////////////////////
+	/** @name Print and IO methods
+	 */
+	//@{
+public:
+	std::string to_string() const {
+		std::ostringstream stream;
+		stream << num_vertices() << " vertices:\n" << vertices_to_string() << std::endl;
+		stream << num_edges() << " edges:\n" << edges_to_string() << std::endl;
+		stream << num_blockers() << " blockers:\n" << blockers_to_string() << std::endl;
+		return stream.str();
+	}
+
+	std::string vertices_to_string() const {
+		std::ostringstream stream;
+		for (auto vertex : vertex_range()) {
+			stream << "{" << (*this)[vertex].get_id() << "} ";
+		}
+		stream << std::endl;
+		return stream.str();
+	}
+
+	std::string edges_to_string() const {
+		std::ostringstream stream;
+		for (auto edge : edge_range())
+			stream << "{" << (*this)[edge].first() << "," << (*this)[edge].second() << "} ";
+		stream << std::endl;
+		return stream.str();
+	}
+
+	std::string blockers_to_string() const {
+		std::ostringstream stream;
+
+		for (auto b : const_blocker_range())
+			stream << *b << std::endl;
+		return stream.str();
+	}
+	//@}
 
 
 
@@ -1533,78 +1544,13 @@ public:
  */
 template<typename Complex, typename SimplexHandleIterator>
 Complex make_complex_from_top_faces(SimplexHandleIterator simplex_begin, SimplexHandleIterator simplex_end, bool is_flag_complex = false) {
-  typedef typename Complex::Simplex_handle Simplex_handle;
-  typedef typename Complex::Blocker_handle Blocker_handle;
-  typedef typename Complex::Vertex_handle Vertex_handle;
-  Complex complex;
-
-  complex.clear();
-
-  std::vector<std::pair<Vertex_handle, Vertex_handle>> edges;
-  // first pass to add vertices and edges
-  for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) {
-    // if meet simplex 9 12 15, need to have at least 16 vertices
-    int max_vertex = 0;
-    for (auto vh : *s_it)
-      max_vertex = (std::max)(max_vertex, static_cast<int>(vh));
-    while (complex.num_vertices() <= max_vertex)
-      complex.add_vertex();
-
-    // for all pairs in s, add an edge
-    for (auto u_it = s_it->begin(); u_it != s_it->end(); ++u_it)
-      for (auto v_it = u_it; ++v_it != s_it->end(); /**/)
-        complex.add_edge(*u_it, *v_it);
-  }
-
-  if (!is_flag_complex) {
-    // need to compute blockers
-
-    // store a structure to decide faster if a simplex is in the complex defined by the max faces
-
-    std::vector<std::set < Simplex_handle >> vertex_to_maxfaces(complex.num_vertices());
-    for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it)
-      vertex_to_maxfaces[s_it->first_vertex()].insert(*s_it);
-
-    // for every maximal face s, it picks its first vertex v and check for all nv \in N(v)
-    // if s union nv is in the complex, if not it is a blocker.
-    for (auto max_face = simplex_begin; max_face != simplex_end; ++max_face) {
-      if (max_face->dimension() > 0) {
-        auto first_v = max_face->first_vertex();
-        for (auto nv : complex.vertex_range(first_v)) {
-          if (!max_face->contains(nv) && max_face->first_vertex() < nv) {
-            // check that all edges in vertices(max_face)\cup nv are here
-            // since max_face is a simplex, we only need to check that edges
-            // (x nv) where x \in max_face are present
-            bool all_edges_here = true;
-            for (auto x : *max_face)
-              if (!complex.contains_edge(x, nv)) {
-                all_edges_here = false;
-                break;
-              }
-            if (all_edges_here) {  // eg  this->contains(max_face)
-              max_face->add_vertex(nv);
-              if (vertex_to_maxfaces[first_v].find(*max_face) == vertex_to_maxfaces[first_v].end()) {  // xxxx
-                // if there exists a blocker included in max_face, we remove it
-                // as it is not a minimum missing face
-                // the other alternative would be to check to all properfaces
-                // are in the complex before adding a blocker but that
-                // would be more expensive if there are few blockers
-                std::vector<Blocker_handle> blockers_to_remove;
-                for (auto b : complex.blocker_range(first_v))
-                  if (b->contains(*max_face))
-                    blockers_to_remove.push_back(b);
-                for (auto b : blockers_to_remove)
-                  complex.delete_blocker(b);
-                complex.add_blocker(*max_face);
-              }
-              max_face->remove_vertex(nv);
-            }
-          }
-        }
-      }
-    }
-  }
-  return complex;
+	typedef typename Complex::Simplex_handle Simplex_handle;
+	std::vector<Simplex_handle> simplices;
+	for (auto top_face = simplex_begin; top_face != simplex_end; ++top_face) {
+		auto subfaces_topface = subfaces(*top_face);
+		simplices.insert(simplices.end(),subfaces_topface.begin(),subfaces_topface.end());
+	}
+	return Complex(simplices.begin(),simplices.end(),is_flag_complex);
 }
 
 }  // namespace skbl
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
index f2c1e447..bb58d0dc 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
@@ -37,98 +37,162 @@ namespace skbl {
  */
 template<typename SkeletonBlockerGeometricDS>
 class Skeleton_blocker_geometric_complex :
-public Skeleton_blocker_complex<SkeletonBlockerGeometricDS> {
- public:
-  typedef typename SkeletonBlockerGeometricDS::GT GT;
-
-  typedef Skeleton_blocker_complex<SkeletonBlockerGeometricDS> SimplifiableSkeletonblocker;
-
-  typedef typename SimplifiableSkeletonblocker::Vertex_handle Vertex_handle;
-  typedef typename SimplifiableSkeletonblocker::Root_vertex_handle Root_vertex_handle;
-  typedef typename SimplifiableSkeletonblocker::Edge_handle Edge_handle;
-  typedef typename SimplifiableSkeletonblocker::Simplex_handle Simplex_handle;
-
-  typedef typename SimplifiableSkeletonblocker::Graph_vertex Graph_vertex;
-
-  typedef typename SkeletonBlockerGeometricDS::Point Point;
-
-  /**
-   * @brief Add a vertex to the complex with its associated point.
-   */
-  Vertex_handle add_vertex(const Point& point) {
-    Vertex_handle ad = SimplifiableSkeletonblocker::add_vertex();
-    (*this)[ad].point() = point;
-    return ad;
-  }
-
-  /**
-   * @brief Returns the Point associated to the vertex v.
-   */
-  const Point& point(Vertex_handle v) const {
-    assert(this->contains_vertex(v));
-    return (*this)[v].point();
-  }
-
-  /**
-   * @brief Returns the Point associated to the vertex v.
-   */
-  Point& point(Vertex_handle v) {
-    assert(this->contains_vertex(v));
-    return (*this)[v].point();
-  }
-
-  const Point& point(Root_vertex_handle global_v) const {
-    Vertex_handle local_v((*this)[global_v]);
-    assert(this->contains_vertex(local_v));
-    return (*this)[local_v].point();
-  }
-
-  Point& point(Root_vertex_handle global_v) {
-    Vertex_handle local_v((*this)[global_v]);
-    assert(this->contains_vertex(local_v));
-    return (*this)[local_v].point();
-  }
-
-  typedef Skeleton_blocker_link_complex<Skeleton_blocker_geometric_complex> Geometric_link;
-
-  /**
-   * Constructs the link of 'simplex' with points coordinates.
-   */
-  Geometric_link link(Vertex_handle v) const {
-    Geometric_link link(*this, Simplex_handle(v));
-    // we now add the point info
-    add_points_to_link(link);
-    return link;
-  }
-
-  /**
-   * Constructs the link of 'simplex' with points coordinates.
-   */
-  Geometric_link link(const Simplex_handle& simplex) const {
-    Geometric_link link(*this, simplex);
-    // we now add the point info
-    add_points_to_link(link);
-    return link;
-  }
-
-  /**
-   * Constructs the link of 'simplex' with points coordinates.
-   */
-  Geometric_link link(Edge_handle edge) const {
-    Geometric_link link(*this, edge);
-    // we now add the point info
-    add_points_to_link(link);
-    return link;
-  }
-
- private:
-
-  void add_points_to_link(Geometric_link& link) const {
-    for (Vertex_handle v : link.vertex_range()) {
-      Root_vertex_handle v_root(link.get_id(v));
-      link.point(v) = (*this).point(v_root);
-    }
-  }
+		public Skeleton_blocker_complex<SkeletonBlockerGeometricDS> {
+		public:
+	typedef typename SkeletonBlockerGeometricDS::GT GT;
+
+	typedef Skeleton_blocker_complex<SkeletonBlockerGeometricDS> SimplifiableSkeletonblocker;
+
+	typedef typename SimplifiableSkeletonblocker::Vertex_handle Vertex_handle;
+	typedef typename SimplifiableSkeletonblocker::Root_vertex_handle Root_vertex_handle;
+	typedef typename SimplifiableSkeletonblocker::Edge_handle Edge_handle;
+	typedef typename SimplifiableSkeletonblocker::Simplex_handle Simplex_handle;
+
+	typedef typename SimplifiableSkeletonblocker::Graph_vertex Graph_vertex;
+
+	typedef typename SkeletonBlockerGeometricDS::Point Point;
+
+
+
+	Skeleton_blocker_geometric_complex(){
+	}
+
+
+
+	/**
+	 * constructor given a list of points
+	 */
+	template<typename PointIterator>
+	explicit Skeleton_blocker_geometric_complex(int num_vertices,PointIterator begin,PointIterator end){
+		for(auto point = begin; point != end; ++point)
+			add_vertex(*point);
+	}
+
+
+	/**
+	 * @brief Constructor with a list of simplices.
+	 * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not).
+	 */
+	template<typename SimpleHandleOutputIterator,typename PointIterator>
+	Skeleton_blocker_geometric_complex(
+			SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end,
+			PointIterator points_begin,PointIterator points_end,
+			bool is_flag_complex = false):Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin,simplex_end,is_flag_complex){
+		unsigned current = 0;
+		for(auto point = points_begin; point != points_end; ++point)
+			(*this)[Vertex_handle(current++)].point() = Point(point->begin(),point->end());
+	}
+
+	template<typename SimpleHandleOutputIterator,typename PointIterator>
+	friend Skeleton_blocker_geometric_complex make_complex_from_top_faces(
+			SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end,
+			PointIterator points_begin,PointIterator points_end,
+			bool is_flag_complex = false){
+		Skeleton_blocker_geometric_complex complex;
+		unsigned current = 0;
+		complex=make_complex_from_top_faces<Skeleton_blocker_geometric_complex>(simplex_begin,simplex_end,is_flag_complex);
+		for(auto point = points_begin; point != points_end; ++point)
+			complex.point(Vertex_handle(current++)) = Point(point->begin(),point->end());
+		return complex;
+	}
+
+
+	/**
+	 * @brief Constructor with a list of simplices.
+	 * Points of every vertex are the point constructed with default constructor.
+	 * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not).
+	 */
+	template<typename SimpleHandleOutputIterator>
+	Skeleton_blocker_geometric_complex(
+			SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end,
+			bool is_flag_complex = false):Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin,simplex_end,is_flag_complex){
+	}
+
+
+	/**
+	 * @brief Add a vertex to the complex with a default constructed associated point.
+	 */
+	Vertex_handle add_vertex() {
+		return SimplifiableSkeletonblocker::add_vertex();
+	}
+
+	/**
+	 * @brief Add a vertex to the complex with its associated point.
+	 */
+	Vertex_handle add_vertex(const Point& point) {
+		Vertex_handle ad = SimplifiableSkeletonblocker::add_vertex();
+		(*this)[ad].point() = point;
+		return ad;
+	}
+
+	/**
+	 * @brief Returns the Point associated to the vertex v.
+	 */
+	const Point& point(Vertex_handle v) const {
+		assert(this->contains_vertex(v));
+		return (*this)[v].point();
+	}
+
+	/**
+	 * @brief Returns the Point associated to the vertex v.
+	 */
+	Point& point(Vertex_handle v) {
+		assert(this->contains_vertex(v));
+		return (*this)[v].point();
+	}
+
+	const Point& point(Root_vertex_handle global_v) const {
+		Vertex_handle local_v((*this)[global_v]);
+		assert(this->contains_vertex(local_v));
+		return (*this)[local_v].point();
+	}
+
+	Point& point(Root_vertex_handle global_v) {
+		Vertex_handle local_v((*this)[global_v]);
+		assert(this->contains_vertex(local_v));
+		return (*this)[local_v].point();
+	}
+
+	typedef Skeleton_blocker_link_complex<Skeleton_blocker_geometric_complex> Geometric_link;
+
+	/**
+	 * Constructs the link of 'simplex' with points coordinates.
+	 */
+	Geometric_link link(Vertex_handle v) const {
+		Geometric_link link(*this, Simplex_handle(v));
+		// we now add the point info
+		add_points_to_link(link);
+		return link;
+	}
+
+	/**
+	 * Constructs the link of 'simplex' with points coordinates.
+	 */
+	Geometric_link link(const Simplex_handle& simplex) const {
+		Geometric_link link(*this, simplex);
+		// we now add the point info
+		add_points_to_link(link);
+		return link;
+	}
+
+	/**
+	 * Constructs the link of 'simplex' with points coordinates.
+	 */
+	Geometric_link link(Edge_handle edge) const {
+		Geometric_link link(*this, edge);
+		// we now add the point info
+		add_points_to_link(link);
+		return link;
+	}
+
+		private:
+
+	void add_points_to_link(Geometric_link& link) const {
+		for (Vertex_handle v : link.vertex_range()) {
+			Root_vertex_handle v_root(link.get_id(v));
+			link.point(v) = (*this).point(v_root);
+		}
+	}
 };
 
 } // namespace skbl
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 163fb7e3..9eab7f1e 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
@@ -468,6 +468,17 @@ template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::add_simplex(const Simplex_handle& sigma) {
 	assert(!this->contains(sigma));
 	assert(sigma.dimension() > 1);
+
+	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);
 }