Fix of cpplint defaults

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

Former-commit-id: 08c5b926a1cde6c41d59f9c143a74dbbd1b1fb5f

4 files changed, 2051 insertions, 2316 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 c98b0b45..aaaab8b0 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
@@ -24,6 +24,7 @@
 
 #include <string>
 #include <vector>
+#include <map>
 
 #include "gudhi/Off_reader.h"
 
@@ -36,38 +37,36 @@ namespace skbl {
  */
 template<typename Complex>
 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_;
-
-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 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 done() {}
+  Complex& complex_;
+  typedef typename Complex::Vertex_handle Vertex_handle;
+  typedef typename Complex::Point Point;
+
+  const bool 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 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 done() { }
 };
 
 /**
@@ -75,137 +74,127 @@ public:
  */
 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()
-		);
-	}
+  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,bool is_flag = false):valid_(false) {
-		std::ifstream stream(name_file);
-		if (stream.is_open()) {
-			if(is_flag || read_only_points){
-				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_;
+ 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 || read_only_points) {
+        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 x : pt)
-					stream<<x<<" ";
-				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;
-	}
-
+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 x : pt)
+          stream << x << " ";
+        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
 
-
 }  // namespace Gudhi
 
-
 #endif  // SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_SKELETON_BLOCKER_OFF_IO_H_
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
index 5b2c1f27..700830f2 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
@@ -63,1479 +63,1466 @@ namespace skbl {
  */
 template<class SkeletonBlockerDS>
 class Skeleton_blocker_complex {
-	template<class ComplexType> friend class Vertex_iterator;
-	template<class ComplexType> friend class Neighbors_vertices_iterator;
-	template<class ComplexType> friend class Edge_iterator;
-	template<class ComplexType> friend class 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;
-	typedef Trie<Simplex_handle> STrie;
-
-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_) {
-		add_vertex_and_edges(simplex_begin,simplex_end);
-
-		if (!is_flag_complex)
-			// need to compute blockers
-			add_blockers(simplex_begin,simplex_end);
-	}
-
-private:
-	template<typename SimpleHandleOutputIterator>
-	void add_vertex_and_edges(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end){
-		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);
-	}
-
-	template<typename SimpleHandleOutputIterator>
-	void add_blockers(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end){
-		Tries<Simplex_handle> tries(num_vertices(),simplex_begin,simplex_end);
-		tries.init_next_dimension();
-		auto simplices(tries.next_dimension_simplices());
-
-		while(!simplices.empty()){
-			simplices = tries.next_dimension_simplices();
-			for(auto& sigma : simplices){
-				//common_positive_neighbors is the set of vertices u such that
-				//for all s in sigma, us is an edge and u>s
-				Simplex_handle common_positive_neighbors(tries.positive_neighbors(sigma.last_vertex()));
-				for(auto sigma_it = sigma.rbegin(); sigma_it != sigma.rend(); ++sigma_it/**/)
-					if(sigma_it != sigma.rbegin())
-						common_positive_neighbors.intersection(tries.positive_neighbors(*sigma_it));
-
-				for(auto x : common_positive_neighbors){
-					//first test that all edges sx are here for all s in sigma
-					bool all_edges_here = true;
-					for(auto s : sigma)
-						if(!contains_edge(x,s)){
-							all_edges_here = false;
-							break;
-						}
-					if(!all_edges_here) continue;
-
-					// all edges of sigma \cup x are here
-					//we have a blocker if all proper faces of sigma \cup x
-					// are in the complex and if sigma \cup x is not in the complex
-					//the first is equivalent at checking if blocks(sigma \cup x) is true
-					//as all blockers of lower dimension have already been computed
-					sigma.add_vertex(x);
-					if(!tries.contains(sigma)&&!blocks(sigma))
-						add_blocker(sigma);
-					sigma.remove_vertex(x);
-				}
-
-			}
-		}
-	}
-
-
-public:
-
-
-	// 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);
-			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 s : sigma)
-			for (auto blocker : const_blocker_range(s))
-				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
-	 */
-	//@{
-
-	/**
-	 * Returns true iff the blocker 'sigma' is popable.
-	 * To define popable, let us call 'L' the complex that
-	 * consists in the current complex without the blocker 'sigma'.
-	 * A blocker 'sigma' is then "popable" if the link of 'sigma'
-	 * in L is reducible.
-	 *
-	 */
-	bool is_popable_blocker(Blocker_handle sigma) const;
-
-	/**
-	 * Removes all the popable blockers of the complex and delete them.
-	 * @returns the number of popable blockers deleted
-	 */
-	void remove_popable_blockers();
-
-	/**
-	 * Removes all the popable blockers of the complex passing through v and delete them.
-	 */
-	void remove_popable_blockers(Vertex_handle v);
-
-	/**
-	 * @brief Removes all the popable blockers of the complex passing through v and delete them.
-	 * Also remove popable blockers in the neighborhood if they became popable.
-	 *
-	 */
-	void remove_all_popable_blockers(Vertex_handle v);
-
-	/**
-	 * Remove the star of the vertex 'v'
-	 */
-	void remove_star(Vertex_handle v) ;
-
-private:
-	/**
-	 * after removing the star of a simplex, blockers sigma that contains this simplex must be removed.
-	 * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added
-	 * whenever the dimension of tau is at least 2.
-	 */
-	void update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed);
-
-public:
-	/**
-	 * Remove the star of the edge connecting vertices a and b.
-	 * @returns the number of blocker that have been removed
-	 */
-	void remove_star(Vertex_handle a, Vertex_handle b) ;
-
-	/**
-	 * Remove the star of the edge 'e'.
-	 */
-	void remove_star(Edge_handle e) ;
-
-	/**
-	 * Remove the star of the simplex 'sigma' which needs to belong to the complex
-	 */
-	void remove_star(const Simplex_handle& sigma);
-
-	/**
-	 * @brief add a maximal simplex plus all its cofaces.
-	 * @details the simplex must have dimension greater than one (otherwise use add_vertex or add_edge).
-	 */
-	void add_simplex(const Simplex_handle& sigma);
-
-private:
-	/**
-	 * remove all blockers that contains sigma
-	 */
-	void remove_blocker_containing_simplex(const Simplex_handle& sigma) ;
-
-	/**
-	 * remove all blockers that contains sigma
-	 */
-	void remove_blocker_include_in_simplex(const Simplex_handle& sigma);
-
-public:
-	enum simplifiable_status {
-		NOT_HOMOTOPY_EQ, MAYBE_HOMOTOPY_EQ, HOMOTOPY_EQ
-	};
-
-	simplifiable_status is_remove_star_homotopy_preserving(const Simplex_handle& simplex) {
-		// todo write a virtual method 'link' in Skeleton_blocker_complex which will be overloaded by the current one of Skeleton_blocker_geometric_complex
-		// then call it there to build the link and return the value of link.is_contractible()
-		return MAYBE_HOMOTOPY_EQ;
-	}
-
-	enum contractible_status {
-		NOT_CONTRACTIBLE, MAYBE_CONTRACTIBLE, CONTRACTIBLE
-	};
-
-	/**
-	 * @brief %Test if the complex is reducible using a strategy defined in the class
-	 * (by default it tests if the complex is a cone)
-	 * @details Note that NO could be returned if some invariant ensures that the complex
-	 * is not a point (for instance if the euler characteristic is different from 1).
-	 * This function will surely have to return MAYBE in some case because the
-	 * associated problem is undecidable but it in practice, it can often
-	 * be solved with the help of geometry.
-	 */
-	virtual contractible_status is_contractible() const {
-		if (this->is_cone()) {
-			return CONTRACTIBLE;
-		} else {
-			return MAYBE_CONTRACTIBLE;
-		}
-	}
-	//@}
-
-	/** @Edge contraction operations
-	 */
-	//@{
-
-	/**
-	 * @return If ignore_popable_blockers is true
-	 * then the result is true iff the link condition at edge ab is satisfied
-	 * or equivalently iff no blocker contains ab.
-	 * If ignore_popable_blockers is false then the
-	 * result is true iff all blocker containing ab are popable.
-	 */
-	bool link_condition(Vertex_handle a, Vertex_handle b, bool ignore_popable_blockers = false) const{
-		for (auto blocker : this->const_blocker_range(a))
-			if (blocker->contains(b)) {
-				// false if ignore_popable_blockers is false
-				// otherwise the blocker has to be popable
-				return ignore_popable_blockers && is_popable_blocker(blocker);
-			}
-		return true;
-
-	}
-
-	/**
-	 * @return If ignore_popable_blockers is true
-	 * then the result is true iff the link condition at edge ab is satisfied
-	 * or equivalently iff no blocker contains ab.
-	 * If ignore_popable_blockers is false then the
-	 * result is true iff all blocker containing ab are popable.
-	 */
-	bool link_condition(Edge_handle e, bool ignore_popable_blockers = false) const {
-		const Graph_edge& edge = (*this)[e];
-		assert(this->get_address(edge.first()));
-		assert(this->get_address(edge.second()));
-		Vertex_handle a(*this->get_address(edge.first()));
-		Vertex_handle b(*this->get_address(edge.second()));
-		return link_condition(a, b, ignore_popable_blockers);
-	}
-
-protected:
-	/**
-	 * Compute simplices beta such that a.beta is an order 0 blocker
-	 * that may be used to construct a new blocker after contracting ab.
-	 * It requires that the link condition is satisfied.
-	 */
-	void tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const;
-
-private:
-	/**
-	 * @brief "Replace" the edge 'bx' by the edge 'ax'.
-	 * Assume that the edge 'bx' was present whereas 'ax' was not.
-	 * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'.
-	 * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added
-	 * and just before edge 'bx' had been removed. That way, it can
-	 * eventually access to information of 'ax'.
-	 */
-	void swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x);
-
-private:
-	/**
-	 * @brief removes all blockers passing through the edge 'ab'
-	 */
-	void delete_blockers_around_vertex(Vertex_handle v);
-
-	/**
-	 * @brief removes all blockers passing through the edge 'ab'
-	 */
-	void delete_blockers_around_edge(Vertex_handle a, Vertex_handle b);
-
-public:
-	/**
-	 * Contracts the edge.
-	 * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
-	 * it removes first all blockers passing through 'ab'
-	 */
-	void contract_edge(Edge_handle edge) {
-		contract_edge(this->first_vertex(edge), this->second_vertex(edge));
-	}
-
-	/**
-	 * Contracts the edge connecting vertices a and b.
-	 * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
-	 * it removes first all blockers passing through 'ab'
-	 */
-	void contract_edge(Vertex_handle a, Vertex_handle b);
-
-private:
-	void get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, std::set<Simplex_handle>& blockers_to_add);
-
-	/**
-	 * delete all blockers that passes through a or b
-	 */
-	void delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b);
-	void update_edges_after_contraction(Vertex_handle a, Vertex_handle b) ;
-
-	void notify_changed_edges(Vertex_handle a) ;
-	//@}
-
-public:
-	/////////////////////////////////////////////////////////////////////////////
-	/** @name Vertex iterators
-	 */
-	//@{
-	typedef Vertex_iterator<Skeleton_blocker_complex> Complex_vertex_iterator;
-
-	//
-	// Range over the vertices of the simplicial complex.
-	// Methods .begin() and .end() return a Complex_vertex_iterator.
-	//
-	typedef boost::iterator_range<Complex_vertex_iterator> 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(this);
-		auto end = Complex_vertex_iterator(this, 0);
-		return Complex_vertex_range(begin, end);
-	}
-
-	typedef Neighbors_vertices_iterator<Skeleton_blocker_complex> Complex_neighbors_vertices_iterator;
-
-
-	typedef boost::iterator_range<Complex_neighbors_vertices_iterator> 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(this, v);
-		auto end = Complex_neighbors_vertices_iterator(this, v, 0);
-		return Complex_neighbors_vertices_range(begin, end);
-	}
-
-	//@}
-
-	/** @name Edge iterators
-	 */
-	//@{
-
-	typedef Edge_iterator<Skeleton_blocker_complex> Complex_edge_iterator;
-
-
-	typedef boost::iterator_range<Complex_edge_iterator> 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(this);
-		auto end = Complex_edge_iterator(this, 0);
-		return Complex_edge_range(begin, end);
-	}
-
-
-	typedef Edge_around_vertex_iterator<Skeleton_blocker_complex> Complex_edge_around_vertex_iterator;
-
-
-	typedef boost::iterator_range <Complex_edge_around_vertex_iterator>	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(this, v);
-		auto end = Complex_edge_around_vertex_iterator(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;
-
-
-	typedef Triangle_iterator<Skeleton_blocker_complex> Complex_triangle_iterator;
-
-
-	/**
-	 * @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();
-	}
-	//@}
-
-
-
-
-
+  template<class ComplexType> friend class Vertex_iterator;
+  template<class ComplexType> friend class Neighbors_vertices_iterator;
+  template<class ComplexType> friend class Edge_iterator;
+  template<class ComplexType> friend class 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;
+  typedef Trie<Simplex_handle> STrie;
+
+ 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_) {
+    add_vertex_and_edges(simplex_begin, simplex_end);
+
+    if (!is_flag_complex)
+      // need to compute blockers
+      add_blockers(simplex_begin, simplex_end);
+  }
+
+ private:
+  template<typename SimpleHandleOutputIterator>
+  void add_vertex_and_edges(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) {
+    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);
+  }
+
+  template<typename SimpleHandleOutputIterator>
+  void add_blockers(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) {
+    Tries<Simplex_handle> tries(num_vertices(), simplex_begin, simplex_end);
+    tries.init_next_dimension();
+    auto simplices(tries.next_dimension_simplices());
+
+    while (!simplices.empty()) {
+      simplices = tries.next_dimension_simplices();
+      for (auto& sigma : simplices) {
+        // common_positive_neighbors is the set of vertices u such that
+        // for all s in sigma, us is an edge and u>s
+        Simplex_handle common_positive_neighbors(tries.positive_neighbors(sigma.last_vertex()));
+        for (auto sigma_it = sigma.rbegin(); sigma_it != sigma.rend(); ++sigma_it)
+          if (sigma_it != sigma.rbegin())
+            common_positive_neighbors.intersection(tries.positive_neighbors(*sigma_it));
+
+        for (auto x : common_positive_neighbors) {
+          // first test that all edges sx are here for all s in sigma
+          bool all_edges_here = true;
+          for (auto s : sigma)
+            if (!contains_edge(x, s)) {
+              all_edges_here = false;
+              break;
+            }
+          if (!all_edges_here) continue;
+
+          // all edges of sigma \cup x are here
+          // we have a blocker if all proper faces of sigma \cup x
+          // are in the complex and if sigma \cup x is not in the complex
+          // the first is equivalent at checking if blocks(sigma \cup x) is true
+          // as all blockers of lower dimension have already been computed
+          sigma.add_vertex(x);
+          if (!tries.contains(sigma) && !blocks(sigma))
+            add_blocker(sigma);
+          sigma.remove_vertex(x);
+        }
+      }
+    }
+  }
+
+ public:
+  // 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);
+      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 s : sigma)
+      for (auto blocker : const_blocker_range(s))
+        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
+   */
+  //@{
+
+  /**
+   * Returns true iff the blocker 'sigma' is popable.
+   * To define popable, let us call 'L' the complex that
+   * consists in the current complex without the blocker 'sigma'.
+   * A blocker 'sigma' is then "popable" if the link of 'sigma'
+   * in L is reducible.
+   *
+   */
+  bool is_popable_blocker(Blocker_handle sigma) const;
+
+  /**
+   * Removes all the popable blockers of the complex and delete them.
+   * @returns the number of popable blockers deleted
+   */
+  void remove_popable_blockers();
+
+  /**
+   * Removes all the popable blockers of the complex passing through v and delete them.
+   */
+  void remove_popable_blockers(Vertex_handle v);
+
+  /**
+   * @brief Removes all the popable blockers of the complex passing through v and delete them.
+   * Also remove popable blockers in the neighborhood if they became popable.
+   *
+   */
+  void remove_all_popable_blockers(Vertex_handle v);
+
+  /**
+   * Remove the star of the vertex 'v'
+   */
+  void remove_star(Vertex_handle v);
+
+ private:
+  /**
+   * after removing the star of a simplex, blockers sigma that contains this simplex must be removed.
+   * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added
+   * whenever the dimension of tau is at least 2.
+   */
+  void update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed);
+
+ public:
+  /**
+   * Remove the star of the edge connecting vertices a and b.
+   * @returns the number of blocker that have been removed
+   */
+  void remove_star(Vertex_handle a, Vertex_handle b);
+
+  /**
+   * Remove the star of the edge 'e'.
+   */
+  void remove_star(Edge_handle e);
+
+  /**
+   * Remove the star of the simplex 'sigma' which needs to belong to the complex
+   */
+  void remove_star(const Simplex_handle& sigma);
+
+  /**
+   * @brief add a maximal simplex plus all its cofaces.
+   * @details the simplex must have dimension greater than one (otherwise use add_vertex or add_edge).
+   */
+  void add_simplex(const Simplex_handle& sigma);
+
+ private:
+  /**
+   * remove all blockers that contains sigma
+   */
+  void remove_blocker_containing_simplex(const Simplex_handle& sigma);
+
+  /**
+   * remove all blockers that contains sigma
+   */
+  void remove_blocker_include_in_simplex(const Simplex_handle& sigma);
+
+ public:
+  enum simplifiable_status {
+    NOT_HOMOTOPY_EQ, MAYBE_HOMOTOPY_EQ, HOMOTOPY_EQ
+  };
+
+  simplifiable_status is_remove_star_homotopy_preserving(const Simplex_handle& simplex) {
+    // todo write a virtual method 'link' in Skeleton_blocker_complex which will be overloaded by the current one of
+    // Skeleton_blocker_geometric_complex
+    // then call it there to build the link and return the value of link.is_contractible()
+    return MAYBE_HOMOTOPY_EQ;
+  }
+
+  enum contractible_status {
+    NOT_CONTRACTIBLE, MAYBE_CONTRACTIBLE, CONTRACTIBLE
+  };
+
+  /**
+   * @brief %Test if the complex is reducible using a strategy defined in the class
+   * (by default it tests if the complex is a cone)
+   * @details Note that NO could be returned if some invariant ensures that the complex
+   * is not a point (for instance if the euler characteristic is different from 1).
+   * This function will surely have to return MAYBE in some case because the
+   * associated problem is undecidable but it in practice, it can often
+   * be solved with the help of geometry.
+   */
+  virtual contractible_status is_contractible() const {
+    if (this->is_cone()) {
+      return CONTRACTIBLE;
+    } else {
+      return MAYBE_CONTRACTIBLE;
+    }
+  }
+  //@}
+
+  /** @Edge contraction operations
+   */
+  //@{
+
+  /**
+   * @return If ignore_popable_blockers is true
+   * then the result is true iff the link condition at edge ab is satisfied
+   * or equivalently iff no blocker contains ab.
+   * If ignore_popable_blockers is false then the
+   * result is true iff all blocker containing ab are popable.
+   */
+  bool link_condition(Vertex_handle a, Vertex_handle b, bool ignore_popable_blockers = false) const {
+    for (auto blocker : this->const_blocker_range(a))
+      if (blocker->contains(b)) {
+        // false if ignore_popable_blockers is false
+        // otherwise the blocker has to be popable
+        return ignore_popable_blockers && is_popable_blocker(blocker);
+      }
+    return true;
+  }
+
+  /**
+   * @return If ignore_popable_blockers is true
+   * then the result is true iff the link condition at edge ab is satisfied
+   * or equivalently iff no blocker contains ab.
+   * If ignore_popable_blockers is false then the
+   * result is true iff all blocker containing ab are popable.
+   */
+  bool link_condition(Edge_handle e, bool ignore_popable_blockers = false) const {
+    const Graph_edge& edge = (*this)[e];
+    assert(this->get_address(edge.first()));
+    assert(this->get_address(edge.second()));
+    Vertex_handle a(*this->get_address(edge.first()));
+    Vertex_handle b(*this->get_address(edge.second()));
+    return link_condition(a, b, ignore_popable_blockers);
+  }
+
+ protected:
+  /**
+   * Compute simplices beta such that a.beta is an order 0 blocker
+   * that may be used to construct a new blocker after contracting ab.
+   * It requires that the link condition is satisfied.
+   */
+  void tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const;
+
+ private:
+  /**
+   * @brief "Replace" the edge 'bx' by the edge 'ax'.
+   * Assume that the edge 'bx' was present whereas 'ax' was not.
+   * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'.
+   * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added
+   * and just before edge 'bx' had been removed. That way, it can
+   * eventually access to information of 'ax'.
+   */
+  void swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x);
+
+ private:
+  /**
+   * @brief removes all blockers passing through the edge 'ab'
+   */
+  void delete_blockers_around_vertex(Vertex_handle v);
+
+  /**
+   * @brief removes all blockers passing through the edge 'ab'
+   */
+  void delete_blockers_around_edge(Vertex_handle a, Vertex_handle b);
+
+ public:
+  /**
+   * Contracts the edge.
+   * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
+   * it removes first all blockers passing through 'ab'
+   */
+  void contract_edge(Edge_handle edge) {
+    contract_edge(this->first_vertex(edge), this->second_vertex(edge));
+  }
+
+  /**
+   * Contracts the edge connecting vertices a and b.
+   * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
+   * it removes first all blockers passing through 'ab'
+   */
+  void contract_edge(Vertex_handle a, Vertex_handle b);
+
+ private:
+  void get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b,
+                                                  std::set<Simplex_handle>& blockers_to_add);
+  /**
+   * delete all blockers that passes through a or b
+   */
+  void delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b);
+  void update_edges_after_contraction(Vertex_handle a, Vertex_handle b);
+  void notify_changed_edges(Vertex_handle a);
+  //@}
+
+ public:
+  /////////////////////////////////////////////////////////////////////////////
+  /** @name Vertex iterators
+   */
+  //@{
+  typedef Vertex_iterator<Skeleton_blocker_complex> Complex_vertex_iterator;
+
+  //
+  // Range over the vertices of the simplicial complex.
+  // Methods .begin() and .end() return a Complex_vertex_iterator.
+  //
+  typedef boost::iterator_range<Complex_vertex_iterator> 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(this);
+    auto end = Complex_vertex_iterator(this, 0);
+    return Complex_vertex_range(begin, end);
+  }
+
+  typedef Neighbors_vertices_iterator<Skeleton_blocker_complex> Complex_neighbors_vertices_iterator;
+
+
+  typedef boost::iterator_range<Complex_neighbors_vertices_iterator> 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(this, v);
+    auto end = Complex_neighbors_vertices_iterator(this, v, 0);
+    return Complex_neighbors_vertices_range(begin, end);
+  }
+
+  //@}
+
+  /** @name Edge iterators
+   */
+  //@{
+
+  typedef Edge_iterator<Skeleton_blocker_complex> Complex_edge_iterator;
+
+
+  typedef boost::iterator_range<Complex_edge_iterator> 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(this);
+    auto end = Complex_edge_iterator(this, 0);
+    return Complex_edge_range(begin, end);
+  }
+
+
+  typedef Edge_around_vertex_iterator<Skeleton_blocker_complex> Complex_edge_around_vertex_iterator;
+
+
+  typedef boost::iterator_range <Complex_edge_around_vertex_iterator> 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(this, v);
+    auto end = Complex_edge_around_vertex_iterator(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;
+  typedef Triangle_iterator<Skeleton_blocker_complex> Complex_triangle_iterator;
+
+  /**
+   * @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();
+  }
+  //@}
 };
 
-
-
 /**
  * build a simplicial complex from a collection
  * of top faces.
  * return the total number of simplices
  */
 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;
-	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);
+Complex make_complex_from_top_faces(SimplexHandleIterator simplex_begin, SimplexHandleIterator simplex_end, 
+                                    bool is_flag_complex = false) {
+  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 69de1085..437482cb 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
@@ -37,195 +37,186 @@ 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;
-
-
-
-	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());
-			complex.point(Vertex_handle(current++)) = Point(*point);
-		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;
-	}
-
-	typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerGeometricDS>> Abstract_link;
-
-	/**
-	 * Constructs the abstract link of v (without points coordinates).
-	 */
-	Abstract_link abstract_link(Vertex_handle v) const {
-		return Abstract_link(*this, Simplex_handle(v));
-	}
-
-	/**
-	 * Constructs the link of 'simplex' with points coordinates.
-	 */
-	Abstract_link abstract_link(const Simplex_handle& simplex) const {
-		return Abstract_link(*this, simplex);
-	}
-
-	/**
-	 * Constructs the link of 'simplex' with points coordinates.
-	 */
-	Abstract_link abstract_link(Edge_handle edge) const {
-		return Abstract_link(*this, edge);
-	}
-
-
-
-		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());
+      complex.point(Vertex_handle(current++)) = Point(*point);
+    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;
+  }
+
+  typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerGeometricDS>> Abstract_link;
+
+  /**
+   * Constructs the abstract link of v (without points coordinates).
+   */
+  Abstract_link abstract_link(Vertex_handle v) const {
+    return Abstract_link(*this, Simplex_handle(v));
+  }
+
+  /**
+   * Constructs the link of 'simplex' with points coordinates.
+   */
+  Abstract_link abstract_link(const Simplex_handle& simplex) const {
+    return Abstract_link(*this, simplex);
+  }
+
+  /**
+   * Constructs the link of 'simplex' with points coordinates.
+   */
+  Abstract_link abstract_link(Edge_handle edge) const {
+    return Abstract_link(*this, edge);
+  }
+
+ 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
+}  // namespace skbl
 
-} // namespace Gudhi
+}  // namespace Gudhi
 
 #endif  // SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_
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 9eab7f1e..86a12d90 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
@@ -32,316 +32,52 @@ namespace Gudhi {
 
 namespace skbl {
 
-///**
-// *  \brief Class that allows simplification operation on a simplicial complex represented
-// *  by a skeleton/blockers pair.
-// * \ingroup skbl
-// * @extends Skeleton_blocker_complex
-// */
-//template<typename SkeletonBlockerDS>
-//class Skeleton_blocker_complex : public Skeleton_blocker_complex<SkeletonBlockerDS> {
-//	template<class ComplexType> friend class Skeleton_blocker_sub_complex;
-//
-//public:
-//	typedef Skeleton_blocker_complex<SkeletonBlockerDS> SkeletonBlockerComplex;
-//
-//	typedef typename SkeletonBlockerComplex::Graph_edge Graph_edge;
-//
-//	typedef typename SkeletonBlockerComplex::boost_adjacency_iterator boost_adjacency_iterator;
-//	typedef typename SkeletonBlockerComplex::Edge_handle Edge_handle;
-//	typedef typename SkeletonBlockerComplex::boost_vertex_handle boost_vertex_handle;
-//	typedef typename SkeletonBlockerComplex::Vertex_handle Vertex_handle;
-//	typedef typename SkeletonBlockerComplex::Root_vertex_handle Root_vertex_handle;
-//	typedef typename SkeletonBlockerComplex::Simplex_handle Simplex_handle;
-//	typedef typename SkeletonBlockerComplex::Root_simplex_handle Root_simplex_handle;
-//	typedef typename SkeletonBlockerComplex::Blocker_handle Blocker_handle;
-//
-//
-//	typedef typename SkeletonBlockerComplex::Root_simplex_iterator Root_simplex_iterator;
-//	typedef typename SkeletonBlockerComplex::Simplex_handle_iterator Simplex_handle_iterator;
-//	typedef typename SkeletonBlockerComplex::BlockerMap BlockerMap;
-//	typedef typename SkeletonBlockerComplex::BlockerPair BlockerPair;
-//	typedef typename SkeletonBlockerComplex::BlockerMapIterator BlockerMapIterator;
-//	typedef typename SkeletonBlockerComplex::BlockerMapConstIterator BlockerMapConstIterator;
-//
-//	typedef typename SkeletonBlockerComplex::Visitor Visitor;
-//
-//
-//	/** @name Constructors / Destructors / Initialization
-//	 */
-//	//@{
-//
-//	explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL) :
-//    				  Skeleton_blocker_complex<SkeletonBlockerDS>(num_vertices_, visitor_) { }
-//
-//	/**
-//	 * @brief Constructor with a list of simplices
-//	 * @details The list of simplices must be the list
-//	 * of simplices of a simplicial complex, sorted with increasing dimension.
-//	 */
-//	//soon deprecated
-//	explicit Skeleton_blocker_complex(std::list<Simplex_handle>& simplices, Visitor* visitor_ = NULL) :
-//    				  Skeleton_blocker_complex<SkeletonBlockerDS>(simplices, visitor_) { }
-//
-//	/**
-//	 * @brief Constructor with a list of simplices
-//	 * @details The list of simplices must be the list of simplices of a simplicial complex.
-//	 */
-//	template<typename SimpleHandleOutputIterator>
-//	explicit Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin,SimpleHandleOutputIterator simplex_end,bool is_flag_complex = false,Visitor* visitor_ = NULL) :
-//	Skeleton_blocker_complex<SkeletonBlockerDS>(simplex_begin,simplex_end,is_flag_complex, visitor_) { }
-//
-//
-//	virtual ~Skeleton_blocker_complex() { }
-//
-//	//@}
-//
-//	/**
-//	 * Returns true iff the blocker 'sigma' is popable.
-//	 * To define popable, let us call 'L' the complex that
-//	 * consists in the current complex without the blocker 'sigma'.
-//	 * A blocker 'sigma' is then "popable" if the link of 'sigma'
-//	 * in L is reducible.
-//	 *
-//	 */
-//	bool is_popable_blocker(Blocker_handle sigma) const;
-//
-//	/**
-//	 * Removes all the popable blockers of the complex and delete them.
-//	 * @returns the number of popable blockers deleted
-//	 */
-//	void remove_popable_blockers();
-//
-//	/**
-//	 * Removes all the popable blockers of the complex passing through v and delete them.
-//	 */
-//	void remove_popable_blockers(Vertex_handle v);
-//
-//	/**
-//	 * @brief Removes all the popable blockers of the complex passing through v and delete them.
-//	 * Also remove popable blockers in the neighborhood if they became popable.
-//	 *
-//	 */
-//	void remove_all_popable_blockers(Vertex_handle v);
-//
-//	/**
-//	 * Remove the star of the vertex 'v'
-//	 */
-//	void remove_star(Vertex_handle v) ;
-//
-//private:
-//	/**
-//	 * after removing the star of a simplex, blockers sigma that contains this simplex must be removed.
-//	 * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added
-//	 * whenever the dimension of tau is at least 2.
-//	 */
-//	void update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed);
-//
-//public:
-//	/**
-//	 * Remove the star of the edge connecting vertices a and b.
-//	 * @returns the number of blocker that have been removed
-//	 */
-//	void remove_star(Vertex_handle a, Vertex_handle b) ;
-//
-//	/**
-//	 * Remove the star of the edge 'e'.
-//	 */
-//	void remove_star(Edge_handle e) ;
-//
-//	/**
-//	 * Remove the star of the simplex 'sigma' which needs to belong to the complex
-//	 */
-//	void remove_star(const Simplex_handle& sigma);
-//
-//	/**
-//	 * @brief add a maximal simplex plus all its cofaces.
-//	 * @details the simplex must have dimension greater than one (otherwise use add_vertex or add_edge).
-//	 */
-//	void add_simplex(const Simplex_handle& sigma);
-//
-//private:
-//	/**
-//	 * remove all blockers that contains sigma
-//	 */
-//	void remove_blocker_containing_simplex(const Simplex_handle& sigma) ;
-//
-//	/**
-//	 * remove all blockers that contains sigma
-//	 */
-//	void remove_blocker_include_in_simplex(const Simplex_handle& sigma);
-//
-//public:
-//	enum simplifiable_status {
-//		NOT_HOMOTOPY_EQ, MAYBE_HOMOTOPY_EQ, HOMOTOPY_EQ
-//	};
-//
-//	simplifiable_status is_remove_star_homotopy_preserving(const Simplex_handle& simplex) {
-//		// todo write a virtual method 'link' in Skeleton_blocker_complex which will be overloaded by the current one of Skeleton_blocker_geometric_complex
-//		// then call it there to build the link and return the value of link.is_contractible()
-//		return MAYBE_HOMOTOPY_EQ;
-//	}
-//
-//	enum contractible_status {
-//		NOT_CONTRACTIBLE, MAYBE_CONTRACTIBLE, CONTRACTIBLE
-//	};
-//
-//	/**
-//	 * @brief %Test if the complex is reducible using a strategy defined in the class
-//	 * (by default it tests if the complex is a cone)
-//	 * @details Note that NO could be returned if some invariant ensures that the complex
-//	 * is not a point (for instance if the euler characteristic is different from 1).
-//	 * This function will surely have to return MAYBE in some case because the
-//	 * associated problem is undecidable but it in practice, it can often
-//	 * be solved with the help of geometry.
-//	 */
-//	virtual contractible_status is_contractible() const {
-//		if (this->is_cone()) {
-//			return CONTRACTIBLE;
-//		} else {
-//			return MAYBE_CONTRACTIBLE;
-//		}
-//	}
-//
-//	/** @Edge contraction operations
-//	 */
-//	//@{
-//
-//	/**
-//	 * @return If ignore_popable_blockers is true
-//	 * then the result is true iff the link condition at edge ab is satisfied
-//	 * or equivalently iff no blocker contains ab.
-//	 * If ignore_popable_blockers is false then the
-//	 * result is true iff all blocker containing ab are popable.
-//	 */
-//	bool link_condition(Vertex_handle a, Vertex_handle b, bool ignore_popable_blockers = false) const{
-//		for (auto blocker : this->const_blocker_range(a))
-//			if (blocker->contains(b)) {
-//				// false if ignore_popable_blockers is false
-//				// otherwise the blocker has to be popable
-//				return ignore_popable_blockers && is_popable_blocker(blocker);
-//			}
-//		return true;
-//
-//	}
-//
-//	/**
-//	 * @return If ignore_popable_blockers is true
-//	 * then the result is true iff the link condition at edge ab is satisfied
-//	 * or equivalently iff no blocker contains ab.
-//	 * If ignore_popable_blockers is false then the
-//	 * result is true iff all blocker containing ab are popable.
-//	 */
-//	bool link_condition(Edge_handle e, bool ignore_popable_blockers = false) const {
-//		const Graph_edge& edge = (*this)[e];
-//		assert(this->get_address(edge.first()));
-//		assert(this->get_address(edge.second()));
-//		Vertex_handle a(*this->get_address(edge.first()));
-//		Vertex_handle b(*this->get_address(edge.second()));
-//		return link_condition(a, b, ignore_popable_blockers);
-//	}
-//
-//protected:
-//	/**
-//	 * Compute simplices beta such that a.beta is an order 0 blocker
-//	 * that may be used to construct a new blocker after contracting ab.
-//	 * It requires that the link condition is satisfied.
-//	 */
-//	void tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const;
-//
-//private:
-//	/**
-//	 * @brief "Replace" the edge 'bx' by the edge 'ax'.
-//	 * Assume that the edge 'bx' was present whereas 'ax' was not.
-//	 * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'.
-//	 * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added
-//	 * and just before edge 'bx' had been removed. That way, it can
-//	 * eventually access to information of 'ax'.
-//	 */
-//	void swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x);
-//
-//private:
-//	/**
-//	 * @brief removes all blockers passing through the edge 'ab'
-//	 */
-//	void delete_blockers_around_vertex(Vertex_handle v);
-//
-//	/**
-//	 * @brief removes all blockers passing through the edge 'ab'
-//	 */
-//	void delete_blockers_around_edge(Vertex_handle a, Vertex_handle b);
-//
-//public:
-//	/**
-//	 * Contracts the edge.
-//	 * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
-//	 * it removes first all blockers passing through 'ab'
-//	 */
-//	void contract_edge(Edge_handle edge) {
-//		contract_edge(this->first_vertex(edge), this->second_vertex(edge));
-//	}
-//
-//	/**
-//	 * Contracts the edge connecting vertices a and b.
-//	 * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
-//	 * it removes first all blockers passing through 'ab'
-//	 */
-//	void contract_edge(Vertex_handle a, Vertex_handle b);
-//
-//private:
-//	void get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, std::set<Simplex_handle>& blockers_to_add);
-//
-//	/**
-//	 * delete all blockers that passes through a or b
-//	 */
-//	void delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b);
-//	void update_edges_after_contraction(Vertex_handle a, Vertex_handle b) ;
-//
-//	void notify_changed_edges(Vertex_handle a) ;
-//	//@}
-//};
-
-
-
+/**
+ * Returns true iff the blocker 'sigma' is popable.
+ * To define popable, let us call 'L' the complex that
+ * consists in the current complex without the blocker 'sigma'.
+ * A blocker 'sigma' is then "popable" if the link of 'sigma'
+ * in L is reducible.
+ *
+ */
 template<typename SkeletonBlockerDS>
 bool Skeleton_blocker_complex<SkeletonBlockerDS>::is_popable_blocker(Blocker_handle sigma) const {
-	assert(this->contains_blocker(*sigma));
-	Skeleton_blocker_link_complex<Skeleton_blocker_complex> link_blocker_sigma;
-	build_link_of_blocker(*this, *sigma, link_blocker_sigma);
+  assert(this->contains_blocker(*sigma));
+  Skeleton_blocker_link_complex<Skeleton_blocker_complex> link_blocker_sigma;
+  build_link_of_blocker(*this, *sigma, link_blocker_sigma);
 
-	bool res = link_blocker_sigma.is_contractible() == CONTRACTIBLE;
-	return res;
+  bool res = link_blocker_sigma.is_contractible() == CONTRACTIBLE;
+  return res;
 }
 
-
 /**
  * Removes all the popable blockers of the complex and delete them.
  * @returns the number of popable blockers deleted
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers() {
-	std::list<Vertex_handle> vertex_to_check;
-	for (auto v : this->vertex_range())
-		vertex_to_check.push_front(v);
-
-	while (!vertex_to_check.empty()) {
-		Vertex_handle v = vertex_to_check.front();
-		vertex_to_check.pop_front();
-
-		bool blocker_popable_found = true;
-		while (blocker_popable_found) {
-			blocker_popable_found = false;
-			for (auto block : this->blocker_range(v)) {
-				if (this->is_popable_blocker(block)) {
-					for (Vertex_handle nv : *block)
-						if (nv != v) vertex_to_check.push_back(nv);
-					this->delete_blocker(block);
-					blocker_popable_found = true;
-					break;
-				}
-			}
-		}
-	}
+  std::list<Vertex_handle> vertex_to_check;
+  for (auto v : this->vertex_range())
+    vertex_to_check.push_front(v);
+
+  while (!vertex_to_check.empty()) {
+    Vertex_handle v = vertex_to_check.front();
+    vertex_to_check.pop_front();
+
+    bool blocker_popable_found = true;
+    while (blocker_popable_found) {
+      blocker_popable_found = false;
+      for (auto block : this->blocker_range(v)) {
+        if (this->is_popable_blocker(block)) {
+          for (Vertex_handle nv : *block)
+            if (nv != v) vertex_to_check.push_back(nv);
+          this->delete_blocker(block);
+          blocker_popable_found = true;
+          break;
+        }
+      }
+    }
+  }
 }
 
 /**
@@ -349,16 +85,16 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers() {
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers(Vertex_handle v) {
-	bool blocker_popable_found = true;
-	while (blocker_popable_found) {
-		blocker_popable_found = false;
-		for (auto block : this->blocker_range(v)) {
-			if (is_popable_blocker(block)) {
-				this->delete_blocker(block);
-				blocker_popable_found = true;
-			}
-		}
-	}
+  bool blocker_popable_found = true;
+  while (blocker_popable_found) {
+    blocker_popable_found = false;
+    for (auto block : this->blocker_range(v)) {
+      if (is_popable_blocker(block)) {
+        this->delete_blocker(block);
+        blocker_popable_found = true;
+      }
+    }
+  }
 }
 
 /**
@@ -368,73 +104,83 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers(Vertex
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_all_popable_blockers(Vertex_handle v) {
-	std::list<Vertex_handle> vertex_to_check;
-	vertex_to_check.push_front(v);
-
-	while (!vertex_to_check.empty()) {
-		Vertex_handle v = vertex_to_check.front();
-		vertex_to_check.pop_front();
-
-		bool blocker_popable_found = true;
-		while (blocker_popable_found) {
-			blocker_popable_found = false;
-			for (auto block : this->blocker_range(v)) {
-				if (this->is_popable_blocker(block)) {
-					for (Vertex_handle nv : *block)
-						if (nv != v) vertex_to_check.push_back(nv);
-					this->delete_blocker(block);
-					blocker_popable_found = true;
-					break;
-				}
-			}
-		}
-	}
+  std::list<Vertex_handle> vertex_to_check;
+  vertex_to_check.push_front(v);
+
+  while (!vertex_to_check.empty()) {
+    Vertex_handle v = vertex_to_check.front();
+    vertex_to_check.pop_front();
+
+    bool blocker_popable_found = true;
+    while (blocker_popable_found) {
+      blocker_popable_found = false;
+      for (auto block : this->blocker_range(v)) {
+        if (this->is_popable_blocker(block)) {
+          for (Vertex_handle nv : *block)
+            if (nv != v) vertex_to_check.push_back(nv);
+          this->delete_blocker(block);
+          blocker_popable_found = true;
+          break;
+        }
+      }
+    }
+  }
 }
 
-
-
+/**
+ * Remove the star of the vertice 'v'
+ */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle v) {
-	// we remove the blockers that are not consistent anymore
-	update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(v));
-	while (this->degree(v) > 0) {
-		Vertex_handle w(* (adjacent_vertices(v.vertex, this->skeleton).first));
-		this->remove_edge(v, w);
-	}
-	this->remove_vertex(v);
+  // we remove the blockers that are not consistent anymore
+  update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(v));
+  while (this->degree(v) > 0) {
+    Vertex_handle w(* (adjacent_vertices(v.vertex, this->skeleton).first));
+    this->remove_edge(v, w);
+  }
+  this->remove_vertex(v);
 }
 
+/**
+ * after removing the star of a simplex, blockers sigma that contains this simplex must be removed.
+ * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added
+ * whenever the dimension of tau is at least 2.
+ */
 template<typename SkeletonBlockerDS>
-void Skeleton_blocker_complex<SkeletonBlockerDS>::update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed) {
-	std::list <Blocker_handle> blockers_to_update;
-	if (simplex_to_be_removed.empty()) return;
-
-	auto v0 = simplex_to_be_removed.first_vertex();
-	for (auto blocker : this->blocker_range(v0)) {
-		if (blocker->contains(simplex_to_be_removed))
-			blockers_to_update.push_back(blocker);
-	}
-
-	for (auto blocker_to_update : blockers_to_update) {
-		Simplex_handle sub_blocker_to_be_added;
-		bool sub_blocker_need_to_be_added =
-				(blocker_to_update->dimension() - simplex_to_be_removed.dimension()) >= 2;
-		if (sub_blocker_need_to_be_added) {
-			sub_blocker_to_be_added = *blocker_to_update;
-			sub_blocker_to_be_added.difference(simplex_to_be_removed);
-		}
-		this->delete_blocker(blocker_to_update);
-		if (sub_blocker_need_to_be_added)
-			this->add_blocker(sub_blocker_to_be_added);
-	}
+void Skeleton_blocker_complex<SkeletonBlockerDS>::update_blockers_after_remove_star_of_vertex_or_edge(
+    const Simplex_handle& simplex_to_be_removed) {
+  std::list <Blocker_handle> blockers_to_update;
+  if (simplex_to_be_removed.empty()) return;
+
+  auto v0 = simplex_to_be_removed.first_vertex();
+  for (auto blocker : this->blocker_range(v0)) {
+    if (blocker->contains(simplex_to_be_removed))
+      blockers_to_update.push_back(blocker);
+  }
+
+  for (auto blocker_to_update : blockers_to_update) {
+    Simplex_handle sub_blocker_to_be_added;
+    bool sub_blocker_need_to_be_added =
+        (blocker_to_update->dimension() - simplex_to_be_removed.dimension()) >= 2;
+    if (sub_blocker_need_to_be_added) {
+      sub_blocker_to_be_added = *blocker_to_update;
+      sub_blocker_to_be_added.difference(simplex_to_be_removed);
+    }
+    this->delete_blocker(blocker_to_update);
+    if (sub_blocker_need_to_be_added)
+      this->add_blocker(sub_blocker_to_be_added);
+  }
 }
 
-
+/**
+ * Remove the star of the edge connecting vertices a and b.
+ * @returns the number of blocker that have been removed
+ */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle a, Vertex_handle b) {
-	update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(a, b));
-	// we remove the edge
-	this->remove_edge(a, b);
+  update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(a, b));
+  // we remove the edge
+  this->remove_edge(a, b);
 }
 
 /**
@@ -442,22 +188,23 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle a, V
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Edge_handle e) {
-	return remove_star(this->first_vertex(e), this->second_vertex(e));
+  return remove_star(this->first_vertex(e), this->second_vertex(e));
 }
 
-
-
+/**
+ * Remove the star of the simplex 'sigma' which needs to belong to the complex
+ */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(const Simplex_handle& sigma) {
-	assert(this->contains(sigma));
-	if (sigma.dimension() == 0) {
-		remove_star(sigma.first_vertex());
-	} else if (sigma.dimension() == 1) {
-		remove_star(sigma.first_vertex(), sigma.last_vertex());
-	} else {
-		remove_blocker_containing_simplex(sigma);
-		this->add_blocker(sigma);
-	}
+  assert(this->contains(sigma));
+  if (sigma.dimension() == 0) {
+    remove_star(sigma.first_vertex());
+  } else if (sigma.dimension() == 1) {
+    remove_star(sigma.first_vertex(), sigma.last_vertex());
+  } else {
+    remove_blocker_containing_simplex(sigma);
+    this->add_blocker(sigma);
+  }
 }
 
 /**
@@ -466,35 +213,34 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(const Simplex_hand
  */
 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);
+  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);
 }
 
-
 /**
  * remove all blockers that contains sigma
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_containing_simplex(const Simplex_handle& sigma) {
-	std::vector <Blocker_handle> blockers_to_remove;
-	for (auto blocker : this->blocker_range(sigma.first_vertex())) {
-		if (blocker->contains(sigma))
-			blockers_to_remove.push_back(blocker);
-	}
-	for (auto blocker_to_update : blockers_to_remove)
-		this->delete_blocker(blocker_to_update);
+  std::vector <Blocker_handle> blockers_to_remove;
+  for (auto blocker : this->blocker_range(sigma.first_vertex())) {
+    if (blocker->contains(sigma))
+      blockers_to_remove.push_back(blocker);
+  }
+  for (auto blocker_to_update : blockers_to_remove)
+    this->delete_blocker(blocker_to_update);
 }
 
 /**
@@ -502,173 +248,195 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_containing_simp
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_include_in_simplex(const Simplex_handle& sigma) {
-	std::vector <Blocker_handle> blockers_to_remove;
-	for (auto blocker : this->blocker_range(sigma.first_vertex())) {
-		if (sigma.contains(*blocker))
-			blockers_to_remove.push_back(blocker);
-	}
-	for (auto blocker_to_update : blockers_to_remove)
-		this->delete_blocker(blocker_to_update);
+  std::vector <Blocker_handle> blockers_to_remove;
+  for (auto blocker : this->blocker_range(sigma.first_vertex())) {
+    if (sigma.contains(*blocker))
+      blockers_to_remove.push_back(blocker);
+  }
+  for (auto blocker_to_update : blockers_to_remove)
+    this->delete_blocker(blocker_to_update);
 }
 
-
+/**
+ * Compute simplices beta such that a.beta is an order 0 blocker
+ * that may be used to construct a new blocker after contracting ab.
+ * It requires that the link condition is satisfied.
+ */
 template<typename SkeletonBlockerDS>
-void Skeleton_blocker_complex<SkeletonBlockerDS>::tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const {
-	for (auto const & blocker : this->const_blocker_range(a)) {
-		Simplex_handle beta = (*blocker);
-		beta.remove_vertex(a);
-		buffer.push_back(beta);
-	}
-
-	Simplex_handle n;
-	this->add_neighbours(b, n);
-	this->remove_neighbours(a, n);
-	n.remove_vertex(a);
-
-
-	for (Vertex_handle y : n) {
-		Simplex_handle beta;
-		beta.add_vertex(y);
-		buffer.push_back(beta);
-	}
+void Skeleton_blocker_complex<SkeletonBlockerDS>::tip_blockers(Vertex_handle a, Vertex_handle b,
+                                                               std::vector<Simplex_handle> & buffer) const {
+  for (auto const & blocker : this->const_blocker_range(a)) {
+    Simplex_handle beta = (*blocker);
+    beta.remove_vertex(a);
+    buffer.push_back(beta);
+  }
+
+  Simplex_handle n;
+  this->add_neighbours(b, n);
+  this->remove_neighbours(a, n);
+  n.remove_vertex(a);
+
+
+  for (Vertex_handle y : n) {
+    Simplex_handle beta;
+    beta.add_vertex(y);
+    buffer.push_back(beta);
+  }
 }
 
+/**
+ * @brief "Replace" the edge 'bx' by the edge 'ax'.
+ * Assume that the edge 'bx' was present whereas 'ax' was not.
+ * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'.
+ * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added
+ * and just before edge 'bx' had been removed. That way, it can
+ * eventually access to information of 'ax'.
+ */
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x) {
-	this->add_edge(a, x);
-	if (this->visitor) this->visitor->on_swaped_edge(a, b, x);
-	this->remove_edge(b, x);
+  this->add_edge(a, x);
+  if (this->visitor) this->visitor->on_swaped_edge(a, b, x);
+  this->remove_edge(b, x);
 }
 
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_vertex(Vertex_handle v) {
-	std::list <Blocker_handle> blockers_to_delete;
-	for (auto blocker : this->blocker_range(v)) {
-		blockers_to_delete.push_back(blocker);
-	}
-	while (!blockers_to_delete.empty()) {
-		this->remove_blocker(blockers_to_delete.back());
-		blockers_to_delete.pop_back();
-	}
+  std::list <Blocker_handle> blockers_to_delete;
+  for (auto blocker : this->blocker_range(v)) {
+    blockers_to_delete.push_back(blocker);
+  }
+  while (!blockers_to_delete.empty()) {
+    this->remove_blocker(blockers_to_delete.back());
+    blockers_to_delete.pop_back();
+  }
 }
 
+/**
+ * @brief removes all blockers passing through the edge 'ab'
+ */
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_edge(Vertex_handle a, Vertex_handle b) {
-	std::list<Blocker_handle> blocker_to_delete;
-	for (auto blocker : this->blocker_range(a))
-		if (blocker->contains(b)) blocker_to_delete.push_back(blocker);
-	while (!blocker_to_delete.empty()) {
-		this->delete_blocker(blocker_to_delete.back());
-		blocker_to_delete.pop_back();
-	}
+  std::list<Blocker_handle> blocker_to_delete;
+  for (auto blocker : this->blocker_range(a))
+    if (blocker->contains(b)) blocker_to_delete.push_back(blocker);
+  while (!blocker_to_delete.empty()) {
+    this->delete_blocker(blocker_to_delete.back());
+    blocker_to_delete.pop_back();
+  }
 }
 
+/**
+ * Contracts the edge connecting vertices a and b.
+ * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
+ * it removes first all blockers passing through 'ab'
+ */
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::contract_edge(Vertex_handle a, Vertex_handle b) {
-	assert(this->contains_vertex(a));
-	assert(this->contains_vertex(b));
-	assert(this->contains_edge(a, b));
+  assert(this->contains_vertex(a));
+  assert(this->contains_vertex(b));
+  assert(this->contains_edge(a, b));
 
-	// if some blockers passes through 'ab', we remove them.
-	if (!link_condition(a, b))
-		delete_blockers_around_edge(a, b);
+  // if some blockers passes through 'ab', we remove them.
+  if (!link_condition(a, b))
+    delete_blockers_around_edge(a, b);
 
-	std::set<Simplex_handle> blockers_to_add;
+  std::set<Simplex_handle> blockers_to_add;
 
-	get_blockers_to_be_added_after_contraction(a, b, blockers_to_add);
+  get_blockers_to_be_added_after_contraction(a, b, blockers_to_add);
 
-	delete_blockers_around_vertices(a, b);
+  delete_blockers_around_vertices(a, b);
 
-	update_edges_after_contraction(a, b);
+  update_edges_after_contraction(a, b);
 
-	this->remove_vertex(b);
+  this->remove_vertex(b);
 
-	notify_changed_edges(a);
+  notify_changed_edges(a);
 
-	for (auto block : blockers_to_add)
-		this->add_blocker(block);
+  for (auto block : blockers_to_add)
+    this->add_blocker(block);
 
-	assert(this->contains_vertex(a));
-	assert(!this->contains_vertex(b));
+  assert(this->contains_vertex(a));
+  assert(!this->contains_vertex(b));
 }
 
 template<typename SkeletonBlockerDS>
 void
-Skeleton_blocker_complex<SkeletonBlockerDS>::get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, std::set<Simplex_handle>& blockers_to_add) {
-	blockers_to_add.clear();
-
-	typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > LinkComplexType;
-
-	LinkComplexType link_a(*this, a);
-	LinkComplexType link_b(*this, b);
-
-	std::vector<Simplex_handle> vector_alpha, vector_beta;
-
-	tip_blockers(a, b, vector_alpha);
-	tip_blockers(b, a, vector_beta);
-
-	for (auto alpha = vector_alpha.begin(); alpha != vector_alpha.end(); ++alpha) {
-		for (auto beta = vector_beta.begin(); beta != vector_beta.end(); ++beta) {
-			Simplex_handle sigma = *alpha;
-			sigma.union_vertices(*beta);
-			Root_simplex_handle sigma_id = this->get_id(sigma);
-			if (this->contains(sigma) &&
-					proper_faces_in_union<Skeleton_blocker_complex<SkeletonBlockerDS>>(sigma_id, link_a, link_b)) {
-				// Blocker_handle blocker = new Simplex_handle(sigma);
-				sigma.add_vertex(a);
-				blockers_to_add.insert(sigma);
-			}
-		}
-	}
+Skeleton_blocker_complex<SkeletonBlockerDS>::get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, 
+    std::set<Simplex_handle>& blockers_to_add) {
+  blockers_to_add.clear();
+
+  typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > LinkComplexType;
+
+  LinkComplexType link_a(*this, a);
+  LinkComplexType link_b(*this, b);
+
+  std::vector<Simplex_handle> vector_alpha, vector_beta;
+
+  tip_blockers(a, b, vector_alpha);
+  tip_blockers(b, a, vector_beta);
+
+  for (auto alpha = vector_alpha.begin(); alpha != vector_alpha.end(); ++alpha) {
+    for (auto beta = vector_beta.begin(); beta != vector_beta.end(); ++beta) {
+      Simplex_handle sigma = *alpha;
+      sigma.union_vertices(*beta);
+      Root_simplex_handle sigma_id = this->get_id(sigma);
+      if (this->contains(sigma) &&
+          proper_faces_in_union<Skeleton_blocker_complex < SkeletonBlockerDS >> (sigma_id, link_a, link_b)) {
+        // Blocker_handle blocker = new Simplex_handle(sigma);
+        sigma.add_vertex(a);
+        blockers_to_add.insert(sigma);
+      }
+    }
+  }
 }
 
+/**
+ * delete all blockers that passes through a or b
+ */
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b) {
-	std::vector<Blocker_handle> blocker_to_delete;
-	for (auto bl : this->blocker_range(a))
-		blocker_to_delete.push_back(bl);
-	for (auto bl : this->blocker_range(b))
-		blocker_to_delete.push_back(bl);
-	while (!blocker_to_delete.empty()) {
-		this->delete_blocker(blocker_to_delete.back());
-		blocker_to_delete.pop_back();
-	}
+  std::vector<Blocker_handle> blocker_to_delete;
+  for (auto bl : this->blocker_range(a))
+    blocker_to_delete.push_back(bl);
+  for (auto bl : this->blocker_range(b))
+    blocker_to_delete.push_back(bl);
+  while (!blocker_to_delete.empty()) {
+    this->delete_blocker(blocker_to_delete.back());
+    blocker_to_delete.pop_back();
+  }
 }
 
-
-
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::update_edges_after_contraction(Vertex_handle a, Vertex_handle b) {
-	// We update the set of edges
-	this->remove_edge(a, b);
-
-	// For all edges {b,x} incident to b,
-	// we remove {b,x} and add {a,x} if not already there.
-	while (this->degree(b) > 0) {
-		Vertex_handle x(*(adjacent_vertices(b.vertex, this->skeleton).first));
-		if (!this->contains_edge(a, x))
-			// we 'replace' the edge 'bx' by the edge 'ax'
-			this->swap_edge(a, b, x);
-		else
-			this->remove_edge(b, x);
-	}
+  // We update the set of edges
+  this->remove_edge(a, b);
+
+  // For all edges {b,x} incident to b,
+  // we remove {b,x} and add {a,x} if not already there.
+  while (this->degree(b) > 0) {
+    Vertex_handle x(*(adjacent_vertices(b.vertex, this->skeleton).first));
+    if (!this->contains_edge(a, x))
+      // we 'replace' the edge 'bx' by the edge 'ax'
+      this->swap_edge(a, b, x);
+    else
+      this->remove_edge(b, x);
+  }
 }
 
 template<typename SkeletonBlockerDS>
 void
-Skeleton_blocker_complex<SkeletonBlockerDS>::notify_changed_edges(Vertex_handle a){
-	// We notify the visitor that all edges incident to 'a' had changed
-	boost_adjacency_iterator v, v_end;
-
-	for (tie(v, v_end) = adjacent_vertices(a.vertex, this->skeleton); v != v_end; ++v)
-		if (this->visitor) this->visitor->on_changed_edge(a, Vertex_handle(*v));
+Skeleton_blocker_complex<SkeletonBlockerDS>::notify_changed_edges(Vertex_handle a) {
+  // We notify the visitor that all edges incident to 'a' had changed
+  boost_adjacency_iterator v, v_end;
 
+  for (tie(v, v_end) = adjacent_vertices(a.vertex, this->skeleton); v != v_end; ++v)
+    if (this->visitor) this->visitor->on_changed_edge(a, Vertex_handle(*v));
 }