From 425b462d361286822ee0ed7b5fe00881ba312ea3 Mon Sep 17 00:00:00 2001
From: vrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb>
Date: Fri, 5 Dec 2014 13:32:54 +0000
Subject: Moved into trunk

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@341 636b058d-ea47-450e-bf9e-a15bfbe3eedb
---
 src/Skeleton_blocker/concept/SkeletonBlockerDS.h   |  127 ++
 .../concept/SkeletonBlockerGeometricDS.h           |   65 +
 src/Skeleton_blocker/doc/SkeletonBlockerMainPage.h |  156 +++
 src/Skeleton_blocker/example/CMakeLists.txt        |    7 +
 .../example/Skeleton_blocker_iteration.cpp         |   73 +
 .../include/gudhi/Skeleton_blocker.h               |   24 +
 .../Skeleton_blocker_complex_visitor.h             |  113 ++
 .../Skeleton_blocker_link_superior.h               |   68 +
 .../Skeleton_blocker/Skeleton_blocker_off_io.h     |  108 ++
 .../Skeleton_blocker_simple_geometric_traits.h     |   66 +
 .../Skeleton_blocker_simple_traits.h               |  142 ++
 .../Skeleton_blocker/Skeleton_blocker_simplex.h    |  375 ++++++
 .../Skeleton_blocker_sub_complex.h                 |  295 ++++
 .../gudhi/Skeleton_blocker/internal/Top_faces.h    |   53 +
 .../Skeleton_blockers_blockers_iterators.h         |  120 ++
 .../iterators/Skeleton_blockers_edges_iterators.h  |  153 +++
 .../iterators/Skeleton_blockers_iterators.h        |   21 +
 .../Skeleton_blockers_simplices_iterators.h        |  412 ++++++
 .../Skeleton_blockers_triangles_iterators.h        |  224 +++
 .../Skeleton_blockers_vertices_iterators.h         |  168 +++
 .../include/gudhi/Skeleton_blocker_complex.h       | 1422 ++++++++++++++++++++
 .../gudhi/Skeleton_blocker_geometric_complex.h     |  118 ++
 .../include/gudhi/Skeleton_blocker_link_complex.h  |  275 ++++
 .../gudhi/Skeleton_blocker_simplifiable_complex.h  |  525 ++++++++
 src/Skeleton_blocker/test/CMakeLists.txt           |    7 +
 src/Skeleton_blocker/test/TestSimplifiable.cpp     |  281 ++++
 .../test/TestSkeletonBlockerComplex.cpp            |  773 +++++++++++
 27 files changed, 6171 insertions(+)
 create mode 100644 src/Skeleton_blocker/concept/SkeletonBlockerDS.h
 create mode 100644 src/Skeleton_blocker/concept/SkeletonBlockerGeometricDS.h
 create mode 100644 src/Skeleton_blocker/doc/SkeletonBlockerMainPage.h
 create mode 100644 src/Skeleton_blocker/example/CMakeLists.txt
 create mode 100644 src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Top_faces.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker_link_complex.h
 create mode 100644 src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
 create mode 100644 src/Skeleton_blocker/test/CMakeLists.txt
 create mode 100644 src/Skeleton_blocker/test/TestSimplifiable.cpp
 create mode 100644 src/Skeleton_blocker/test/TestSkeletonBlockerComplex.cpp

(limited to 'src/Skeleton_blocker')

diff --git a/src/Skeleton_blocker/concept/SkeletonBlockerDS.h b/src/Skeleton_blocker/concept/SkeletonBlockerDS.h
new file mode 100644
index 00000000..1fa9dc06
--- /dev/null
+++ b/src/Skeleton_blocker/concept/SkeletonBlockerDS.h
@@ -0,0 +1,127 @@
+/*
+ * SkeletonBlockerDS.h
+ *
+ *  Created on: Feb 20, 2014 
+ *      Author: David Salinas
+ *  Copyright 2013 INRIA. All rights reserved
+ */
+
+#ifndef GUDHI_SKELETONBLOCKERDS_H_
+#define GUDHI_SKELETONBLOCKERDS_H_
+
+
+namespace Gudhi {
+
+namespace skbl {
+
+
+
+/** \brief Concept that must be passed to
+ * the template class Skeleton_blockers_complex
+ *
+ */
+struct SkeletonBlockerDS
+{
+	/**
+	 * @todo faire un default value pour les vertex_handle
+	 */
+
+	/**
+	 * @brief index that allows to find the vertex in the boost graph
+	 */
+	typedef int boost_vertex_handle;
+
+
+	/**
+	 * @brief Root_vertex_handle and Vertex_handle are similar to global and local vertex descriptor
+	 * used in <a href="http://www.boost.org/doc/libs/1_38_0/libs/graph/doc/subgraph.html">boost subgraphs</a>
+	 * and allow to localize a vertex of a subcomplex on its parent root complex.
+	 *
+	 * In gross, vertices are stored in a vector
+	 * and the Root_vertex_handle and Vertex_handle store indices of a vertex in this vector.
+	 *
+	 * For the root simplicial complex, the Root_vertex_handle and Vertex_handle of a vertex
+	 * are the same.
+	 *
+	 *
+	 * For a subcomplex L of a simplicial complex K, the local descriptor, ie the Vertex_handle, of a
+	 * vertex v (that belongs to L)  is its position in the vector of vertices
+	 * of the subcomplex L whereas its Root_vertex_handle (global descriptor) is the position of v in the vector of the
+	 * vertices of the root simplicial complex K.
+	 */
+	struct Root_vertex_handle{
+
+		boost_vertex_handle vertex;
+
+		friend ostream& operator << (ostream& o, const Root_vertex_handle & v);
+	};
+
+	/**
+	 * A Vertex_handle must be Default Constructible, Assignable and Equality Comparable.
+	 */
+	struct Vertex_handle{
+		boost_vertex_handle vertex;
+
+		friend ostream& operator << (ostream& o, const Vertex_handle & v);
+	};
+
+
+	/**
+	 * \brief The type of vertices that are stored the boost graph.
+	 * A Vertex must be Default Constructible and Equality Comparable.
+	 *
+	 */
+	struct Graph_vertex{
+		/** \brief Used to deactivate a vertex for example when contracting an edge.
+		 * It allows in some cases to remove the vertex at low cost.
+		 */
+		void deactivate();
+
+		/** \brief Used to activate a vertex.
+		 */
+		void activate();
+
+		/** \brief Tells if the vertex is active.
+		 */
+		bool is_active() const;
+
+		void set_id(Root_vertex_handle i);
+		Root_vertex_handle get_id() const;
+		virtual string to_string() const ;
+		friend ostream& operator << (ostream& o, const Graph_vertex & v);
+	};
+
+
+	/**
+	 * \brief The type of edges that are stored the boost graph.
+	 * An Edge must be Default Constructible and Equality Comparable.
+	 */
+	struct Graph_edge{
+		/**
+		 * @brief Allows to modify vertices of the edge.
+		 */
+		void setId(Root_vertex_handle a,Root_vertex_handle b);
+
+		/**
+		 * @brief Returns the first vertex of the edge.
+		 */
+		Root_vertex_handle first() const ;
+
+		/**
+		 * @brief Returns the second vertex of the edge.
+		 */
+		Root_vertex_handle second() const ;
+
+		friend ostream& operator << (ostream& o, const Simple_edge & v);
+	};
+
+
+};
+
+
+
+}  // namespace skbl
+}  // namespace GUDHI
+
+
+#endif /* GUDHI_SKELETONBLOCKERDS_H_ */
diff --git a/src/Skeleton_blocker/concept/SkeletonBlockerGeometricDS.h b/src/Skeleton_blocker/concept/SkeletonBlockerGeometricDS.h
new file mode 100644
index 00000000..23edaaef
--- /dev/null
+++ b/src/Skeleton_blocker/concept/SkeletonBlockerGeometricDS.h
@@ -0,0 +1,65 @@
+/*
+ * SkeletonBlockerGeometricDS.h
+ *
+ *  Created on: Feb 20, 2014 
+ *      Author: David Salinas
+ *  Copyright 2013 INRIA. All rights reserved
+ */
+
+#ifndef GUDHI_SKELETONBLOCKERGEOMETRICDS_H_
+#define GUDHI_SKELETONBLOCKERGEOMETRICDS_H_
+
+/** \brief Concept that must be passed to
+ * the template class Skeleton_blocker_geometric_complex
+ *
+ */
+template<typename GeometryTrait>
+struct SkeletonBlockerGeometricDS : public SkeletonBlockerDS
+{
+
+	/**
+	 * Geometry information.
+	 */
+	typedef GeometryTrait GT ;
+
+	/**
+	 * Type of point (should be the same as GT::Point).
+	 */
+	typedef typename GeometryTrait::Point Point;
+
+	/**
+	 * @brief Vertex that stores a point.
+	 */
+	class Graph_vertex : public SkeletonBlockerDS::Graph_vertex{
+	public:
+		/**
+		 * @brief Access to the point.
+		 */
+		Point& point(){	return point_; }
+		/**
+		 * @brief Access to the point.
+		 */
+		const Point& point() const {	return point_; }
+	};
+
+	/**
+	 * @brief Edge that allows to access to an index.
+	 * The indices of the edges are used to store heap information
+	 * in the edge contraction algorithm.
+	 */
+	class Graph_Edge : public SkeletonBlockerDS::Graph_edge{
+	public:
+		/**
+		 * @brief Access to the index.
+		 */
+		int& index();
+		/**
+		 * @brief Access to the index.
+		 */
+		int index();
+	};
+};
+
+
+
+#endif /* GUDHI_SKELETONBLOCKERGEOMETRICDS_H_ */
diff --git a/src/Skeleton_blocker/doc/SkeletonBlockerMainPage.h b/src/Skeleton_blocker/doc/SkeletonBlockerMainPage.h
new file mode 100644
index 00000000..5a80bc8d
--- /dev/null
+++ b/src/Skeleton_blocker/doc/SkeletonBlockerMainPage.h
@@ -0,0 +1,156 @@
+
+/*! \mainpage Skeleton blockers data-structure
+
+\author David Salinas
+
+\section Introduction
+The Skeleton-blocker data-structure had been introduced in the two papers 
+\cite skbl_socg2011 \cite skbl_ijcga2012. 
+It proposes a light encoding for simplicial complexes by storing only an *implicit* representation of its
+simplices.
+Intuitively, it just stores the 1-skeleton of a simplicial complex with a graph and the set of its "missing faces" that
+is very small in practice (see next section for a formal definition).
+This data-structure handles every classical operations used for simplicial complexes such as
+ as simplex enumeration or simplex removal but operations that are particularly efficient 
+ are operations that do not require simplex enumeration such as edge iteration, link computation or simplex contraction.
+
+
+\todo{image wont include}
+
+\section Definitions
+
+We recall briefly classical definitions of simplicial complexes  \cite Munkres.
+An abstract simplex is a finite non-empty set and its dimension is its number of elements minus 1.
+Whenever \f$\tau \subset \sigma\f$ and \f$\tau \neq \emptyset \f$, \f$ \tau \f$ is called a face of 
+\f$ \sigma\f$  and \f$ \sigma\f$ is called a coface of \f$ \tau \f$ . Furthermore,
+when \f$ \tau \neq \sigma\f$ we say that \f$ \tau\f$ is a proper-face of \f$ \sigma\f$.
+An abstract simplicial complex is a set of simplices that contains all the faces of its simplices.
+The 1-skeleton of a simplicial complex (or its graph) consists of its elements of dimension lower than 2.
+
+
+\image latex "images/ds_representation.eps" "My application" width=10cm
+
+
+To encode, a simplicial complex, one can encodes all its simplices. 
+In case when this number gets too large,
+a lighter and implicit version consists of encoding only its graph plus some elements called missing faces or blockers.
+A blocker is a simplex of dimension greater than 1 
+that does not belong to the complex but whose all proper faces does.
+
+
+Remark that for a clique complex (i.e. a simplicial complex whose simplices are cliques of its graph), the set of blockers
+is empty and the data-structure is then particularly sparse. 
+One famous example of clique-complex is the Rips complex which is intensively used
+in topological data-analysis.
+In practice, the set of blockers of a simplicial complex 
+remains also small when simplifying a Rips complex with edge contractions 
+but also for most of the simplicial complexes used in topological data-analysis such as Delaunay, Cech or Witness complexes. 
+For instance, the numbers of blockers is depicted for a random 3 dimensional sphere embedded into \f$R^4\f$ 
+in figure X.
+
+
+\image latex "images/blockers_curve.eps" "My application" width=10cm
+
+
+
+
+\section API
+
+\subsection Overview
+
+Four classes define  a simplicial complex namely :
+
+\li <Code>Skeleton_blocker_complex</Code> : a simplicial complex with basic operations such as vertex/edge/simplex enumeration and construction
+\li <Code>Skeleton_blocker_link_complex</Code> : the link of a simplex in a parent complex. It is represented as a sub complex
+of the parent complex
+\li <Code>Skeleton_blocker_simplifiable_complex</Code> : a simplicial complex with simplification operations such as edge contraction or simplex collapse
+\li <Code>Skeleton_blocker_geometric_complex</Code> : a simplicial complex who has access to geometric points in  \f$R^d\f$ 
+
+The two last classes are derived classes from the <Code>Skeleton_blocker_complex</Code> class. The class <Code>Skeleton_blocker_link_complex</Code> inheritates from a template passed parameter
+that may be either <Code>Skeleton_blocker_complex</Code> or <Code>Skeleton_blocker_geometric_complex</Code> (a link may store points coordinates or not).
+
+\todo{include links}
+
+\subsection Visitor
+
+The class <Code>Skeleton_blocker_complex</Code> has a visitor that is called when usual operations such as adding an edge or remove a vertex are called.
+You may want to use this visitor to compute statistics or to update another data-structure (for instance this visitor is heavily used in the 
+<Code>Contraction</Code> package).
+
+
+
+\section Example
+
+ 
+\subsection s Iterating through vertices, edges, blockers and simplices	
+ 
+Iteration through vertices, edges, simplices or blockers is straightforward with c++11 for range loops.
+Note that simplex iteration with this implicit data-structure just takes
+a few more time compared to iteration via an explicit representation 
+such as the Simplex Tree. The following example computes the Euler Characteristic
+of a simplicial complex.
+
+  \code{.cpp}
+	typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> Complex;
+	typedef Complex::Vertex_handle Vertex_handle;
+	typedef Complex::Simplex_handle Simplex;
+
+  	const int n = 15;
+
+	// build a full complex with 10 vertices and 2^n-1 simplices
+	Complex complex;
+	for(int i=0;i<n;i++)
+		complex.add_vertex();
+	for(int i=0;i<n;i++)
+		for(int j=0;j<i;j++)
+			//note that add_edge adds the edge and all its cofaces
+			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+
+	// this is just to illustrate iterators, to count number of vertices
+	// or edges, complex.num_vertices() and complex.num_edges() are
+	// more appropriated!
+	unsigned num_vertices = 0;
+	for(auto v : complex.vertex_range()){
+		++num_vertices;
+	}
+
+	unsigned num_edges = 0;
+	for(auto e : complex.edge_range())
+		++num_edges;
+
+	unsigned euler = 0;
+	unsigned num_simplices = 0;
+	// we use a reference to a simplex instead of a copy
+	// value here because a simplex is a set of integers
+	// and copying it cost time
+	for(const Simplex & s : complex.simplex_range()){
+		++num_simplices;
+		if(s.dimension()%2 == 0) 
+			euler += 1;
+		else 
+			euler -= 1;
+	}
+	std::cout << "Saw "<<num_vertices<<" vertices, "<<num_edges<<" edges and "<<num_simplices<<" simplices"<<std::endl;
+	std::cout << "The Euler Characteristic is "<<euler<<std::endl;
+  \endcode
+
+
+\verbatim
+./SkeletonBlockerIteration
+Saw 15 vertices, 105 edges and 32767 simplices
+The Euler Characteristic is 1
+ 0.537302s wall, 0.530000s user + 0.000000s system = 0.530000s CPU (98.6%)
+\endverbatim
+
+
+
+\subsection Acknowledgements
+The author wishes to thank Dominique Attali and André Lieutier for 
+their collaboration to write the two initial papers about this data-structure
+ and also Dominique for leaving him use a prototype. 
+
+
+\copyright GNU General Public License v3.                         
+\verbatim  Contact: David Salinas,     david.salinas@inria.fr \endverbatim
+
+*/
diff --git a/src/Skeleton_blocker/example/CMakeLists.txt b/src/Skeleton_blocker/example/CMakeLists.txt
new file mode 100644
index 00000000..0aed1e7a
--- /dev/null
+++ b/src/Skeleton_blocker/example/CMakeLists.txt
@@ -0,0 +1,7 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIskbl)
+
+add_executable ( SkeletonBlockerIteration Skeleton_blocker_iteration.cpp )
+
+target_link_libraries( SkeletonBlockerIteration ${Boost_TIMER_LIBRARY})	
+
diff --git a/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp b/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp
new file mode 100644
index 00000000..3f4a9ffc
--- /dev/null
+++ b/src/Skeleton_blocker/example/Skeleton_blocker_iteration.cpp
@@ -0,0 +1,73 @@
+#include <boost/timer/timer.hpp>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+#include <fstream>
+#include <sstream>
+
+
+#include "gudhi/Skeleton_blocker.h"
+//#include "gudhi/Skeleton_blocker_complex.h"
+//#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h"
+
+using namespace std;
+using namespace Gudhi;
+using namespace skbl;
+
+typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> Complex;
+typedef Complex::Vertex_handle Vertex_handle;
+typedef Complex::Simplex_handle Simplex;
+
+
+Complex build_complete_complex(int n){
+	// build a full complex with 10 vertices and 2^n-1 simplices
+		Complex complex;
+		for(int i=0;i<n;i++)
+			complex.add_vertex();
+		for(int i=0;i<n;i++)
+			for(int j=0;j<i;j++)
+				//note that add_edge, add the edge and all its cofaces
+				complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+		return complex;
+}
+
+int main (int argc, char *argv[]){
+	boost::timer::auto_cpu_timer t;
+
+	const int n = 15;
+
+	// build a full complex with 10 vertices and 2^n-1 simplices
+	Complex complex(build_complete_complex(n));
+
+	// this is just to illustrate iterators, to count number of vertices
+	// or edges, complex.num_vertices() and complex.num_edges() are
+	// more appropriated!
+	unsigned num_vertices = 0;
+	for(auto v : complex.vertex_range()){
+		if(v==v); //do something with v as removing a ennoying warning.
+		++num_vertices;
+	}
+
+	unsigned num_edges = 0;
+	for(auto e : complex.edge_range()){
+		if(e==e);
+		++num_edges;
+	}
+
+	unsigned euler = 0;
+	unsigned num_simplices = 0;
+	// we use a reference to a simplex instead of a copy
+	// value here because a simplex is a set of integers
+	// and copying it cost time
+	for(const Simplex & s : complex.simplex_range()){
+		++num_simplices;
+		if(s.dimension()%2 == 0) 
+			euler += 1;
+		else 
+			euler -= 1;
+	}
+	std::cout << "Saw "<<num_vertices<<" vertices, "<<num_edges<<" edges and "<<num_simplices<<" simplices"<<std::endl;
+	std::cout << "The Euler Characteristic is "<<euler<<std::endl;
+	return EXIT_SUCCESS;
+}
+
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h
new file mode 100644
index 00000000..1d215984
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h
@@ -0,0 +1,24 @@
+/*
+ * Skeleton_blocker.h
+ *
+ *  Created on: Jan, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKER_H
+#define GUDHI_SKELETON_BLOCKER_H
+
+#include "gudhi/Skeleton_blocker_complex.h"
+#include "gudhi/Skeleton_blocker_geometric_complex.h"
+#include "gudhi/Skeleton_blocker_simplifiable_complex.h"
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h"
+
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h"
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h"
+
+
+#include "gudhi/Utils.h" //xxx
+
+
+#endif
+
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h
new file mode 100644
index 00000000..4475afcc
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h
@@ -0,0 +1,113 @@
+/*
+ * ComplexVisitor.h
+ *
+ *  Created on: Dec 11, 2013
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_COMPLEXVISITOR_H_
+#define GUDHI_COMPLEXVISITOR_H_
+
+
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h"
+
+namespace Gudhi{
+
+namespace skbl {
+// todo rajouter les const
+
+/**
+ *@class Skeleton_blocker_complex_visitor
+ *@brief Interface for a visitor of a simplicial complex.
+ */
+template <typename Vertex_handle>
+class Skeleton_blocker_complex_visitor {
+public:
+	virtual ~Skeleton_blocker_complex_visitor(){};
+
+	virtual void on_add_vertex(Vertex_handle) = 0;
+	virtual void on_remove_vertex(Vertex_handle) = 0;
+
+	virtual void on_add_edge(Vertex_handle a,Vertex_handle b) = 0;
+	virtual void on_remove_edge(Vertex_handle a,Vertex_handle b) = 0;
+
+	/**
+	 * @brief Called when an edge changes, during the contraction of
+	 * an edge
+	 */
+	virtual void on_changed_edge(Vertex_handle a,Vertex_handle b) = 0;
+
+	/**
+	 * @brief Called when performing an edge contraction when
+	 * an edge bx is replaced by an edge ax (not already present).
+	 * Precisely, this methods is called this way in contract_edge :
+	 * add_edge(a,x)
+	 * on_swaped_edge(a,b,x)
+	 * remove_edge(b,x)
+	 */
+	virtual void on_swaped_edge(Vertex_handle a,Vertex_handle b,Vertex_handle x)=0;
+	virtual void on_add_blocker(const Skeleton_blocker_simplex<Vertex_handle>&) = 0;
+	virtual void on_delete_blocker(const Skeleton_blocker_simplex<Vertex_handle>*) = 0;
+};
+
+
+/**
+ *@class Dummy_complex_visitor
+ *@brief A dummy visitor of a simplicial complex that does nothing
+ *
+ */
+template <typename Vertex_handle>
+class Dummy_complex_visitor  : public Skeleton_blocker_complex_visitor<Vertex_handle>{
+public:
+	void on_add_vertex(Vertex_handle) {}
+	void on_remove_vertex(Vertex_handle){}
+	void on_add_edge(Vertex_handle a,Vertex_handle b){}
+	void on_remove_edge(Vertex_handle a,Vertex_handle b){}
+	void on_changed_edge(Vertex_handle a,Vertex_handle b){}
+	void on_swaped_edge(Vertex_handle a,Vertex_handle b,Vertex_handle x){}
+	void on_add_blocker(const Skeleton_blocker_simplex<Vertex_handle>&){}
+	void on_delete_blocker(const Skeleton_blocker_simplex<Vertex_handle>*){}
+
+};
+
+
+
+/**
+ *@class Print_complex_visitor
+ *@brief A visitor that prints the visit information.
+ *
+ */
+template <typename Vertex_handle>
+class Print_complex_visitor  : public  Skeleton_blocker_complex_visitor<Vertex_handle>{
+public:
+	void on_add_vertex(Vertex_handle v) {
+		std::cerr << "on_add_vertex:"<<v<<std::endl;
+	}
+	void on_remove_vertex(Vertex_handle v){
+		std::cerr << "on_remove_vertex:"<<v<<std::endl;
+	}
+	void on_add_edge(Vertex_handle a,Vertex_handle b){
+		std::cerr << "on_add_edge:"<<a<<","<<b<<std::endl;
+	}
+	void on_remove_edge(Vertex_handle a,Vertex_handle b){
+		std::cerr << "on_remove_edge:"<<a<<","<<b<<std::endl;
+	}
+	void on_changed_edge(Vertex_handle a,Vertex_handle b){
+		std::cerr << "on_changed_edge:"<<a<<","<<b<<std::endl;
+	}
+	void on_swaped_edge(Vertex_handle a,Vertex_handle b,Vertex_handle x){
+			std::cerr << "on_swaped_edge:"<<a<<","<<b<<","<<x<<std::endl;
+	}
+	void on_add_blocker(const Skeleton_blocker_simplex<Vertex_handle>& b){
+		std::cerr << "on_add_blocker:"<<b<<std::endl;
+	}
+	void on_delete_blocker(const Skeleton_blocker_simplex<Vertex_handle>* b){
+		std::cerr << "on_delete_blocker:"<<b<<std::endl;
+	}
+};
+
+}
+
+}  // namespace GUDHI
+
+#endif /* GUDHI_COMPLEXVISITOR_H_ */
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h
new file mode 100644
index 00000000..864ee6ac
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h
@@ -0,0 +1,68 @@
+/*
+ * Skeleton_blocker_link_superior.h
+ *
+ *  Created on: Feb 19, 2014 
+ *      Author: David Salinas
+ *  Copyright 2013 INRIA. All rights reserved
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKER_LINK_SUPERIOR_H_
+#define GUDHI_SKELETON_BLOCKER_LINK_SUPERIOR_H_
+
+#include "gudhi/Skeleton_blocker_link_complex.h"
+
+namespace Gudhi{
+
+namespace skbl {
+
+template<class ComplexType> class Skeleton_blocker_sub_complex;
+
+
+/**
+ *  \brief Class representing the link of a simplicial complex encoded by a skeleton/blockers pair.
+ *  It computes only vertices greater than the simplex used to build the link.
+ */
+template<typename ComplexType>
+class Skeleton_blocker_link_superior : public Skeleton_blocker_link_complex<ComplexType>
+{
+	typedef typename ComplexType::Edge_handle Edge_handle;
+
+	typedef typename ComplexType::boost_vertex_handle boost_vertex_handle;
+
+public:
+
+	typedef typename ComplexType::Vertex_handle Vertex_handle;
+	typedef typename ComplexType::Root_vertex_handle Root_vertex_handle;
+	typedef typename ComplexType::Simplex_handle Simplex_handle;
+	typedef typename ComplexType::Root_simplex_handle Root_simplex_handle;
+	typedef typename ComplexType::BlockerMap BlockerMap;
+	typedef typename ComplexType::BlockerPair BlockerPair;
+	typedef typename ComplexType::BlockerMapIterator BlockerMapIterator;
+	typedef typename ComplexType::BlockerMapConstIterator BlockerMapConstIterator;
+	typedef typename ComplexType::Simplex_handle::Simplex_vertex_const_iterator AddressSimplexConstIterator;
+	typedef typename ComplexType::Root_simplex_handle::Simplex_vertex_const_iterator IdSimplexConstIterator;
+
+
+	Skeleton_blocker_link_superior()
+		:Skeleton_blocker_link_complex<ComplexType>(true)
+		 {
+		 }
+
+	Skeleton_blocker_link_superior(const ComplexType & parent_complex, Simplex_handle& alpha_parent_adress)
+		:Skeleton_blocker_link_complex<ComplexType>(parent_complex,alpha_parent_adress,true)
+		 {
+		 }
+
+	Skeleton_blocker_link_superior(const ComplexType & parent_complex, Vertex_handle a_parent_adress)
+	:Skeleton_blocker_link_complex<ComplexType>(parent_complex,a_parent_adress,true)
+	 {
+	}
+
+};
+
+}
+
+}  // namespace GUDHI
+
+
+#endif /* SKELETON_BLOCKER_LINK_SUPERIOR_H_ */
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
new file mode 100644
index 00000000..369635e5
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h
@@ -0,0 +1,108 @@
+/*
+ * Skeleton_blocker_off_io.h
+ *  Created on: Nov 28, 2014
+ * This file is part of the Gudhi Library. The Gudhi library 
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++ 
+ *    library for computational topology.
+ *
+ *    Author(s):       David Salinas
+ *
+ *    Copyright (C) 2014  INRIA Sophia Antipolis-Méditerranée (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ * 
+ */
+
+
+#ifndef SKELETON_BLOCKER_OFF_IO_H_
+#define SKELETON_BLOCKER_OFF_IO_H_
+
+#include "gudhi/Off_reader.h"
+
+namespace Gudhi {
+
+namespace skbl {
+
+template<typename Complex>
+class Skeleton_blocker_off_visitor_reader{
+	Complex& complex_;
+	typedef typename Complex::Vertex_handle Vertex_handle;
+	typedef typename Complex::Point Point;
+
+	const bool load_only_points_;
+
+public:
+	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){
+		//todo do an assert to check that this number are correctly read
+		//todo reserve size for vector points
+	}
+
+
+	void point(const std::vector<double>& point){
+		complex_.add_vertex(point);
+	}
+
+	void 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(){
+	}
+};
+
+template<typename Complex>
+class Skeleton_blocker_off_reader{
+public:
+	/**
+	 * name_file : file to read
+	 * read_complex : complex that will receive the file content
+	 * read_only_points : specify true if only the points must be read
+	 */
+	Skeleton_blocker_off_reader(const std::string & name_file,Complex& read_complex,bool read_only_points = false):valid_(false){
+		std::ifstream stream(name_file);
+		if(stream.is_open()){
+			Skeleton_blocker_off_visitor_reader<Complex> off_visitor(read_complex,read_only_points);
+			Off_reader off_reader(stream);
+			valid_ = off_reader.read(off_visitor);
+		}
+	}
+
+	/**
+	 * return true iff reading did not meet problems.
+	 */
+	bool is_valid() const{
+		return valid_;
+	}
+
+private:
+	bool valid_;
+};
+
+}  // namespace skbl
+
+
+}  // namespace Gudhi
+
+
+#endif /* SKELETON_BLOCKER_OFF_IO_H_ */
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h
new file mode 100644
index 00000000..79a8e7aa
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h
@@ -0,0 +1,66 @@
+/*
+ * Skeleton_blocker_simple_geometric_traits.h
+ *
+ *  Created on: Feb 11, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKERS_SIMPLE_GEOMETRIC_TRAITS_H_
+#define GUDHI_SKELETON_BLOCKERS_SIMPLE_GEOMETRIC_TRAITS_H_
+
+#include <string>
+#include <sstream>
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h"
+
+namespace Gudhi{
+
+namespace skbl{
+
+/**
+ * @extends SkeletonBlockerGeometricDS
+ */
+template<typename GeometryTrait>
+struct Skeleton_blocker_simple_geometric_traits : public skbl::Skeleton_blocker_simple_traits {
+public:
+
+	typedef GeometryTrait GT;
+	typedef typename GT::Point Point;
+	typedef typename Skeleton_blocker_simple_traits::Root_vertex_handle Root_vertex_handle;
+	typedef typename Skeleton_blocker_simple_traits::Graph_vertex Simple_vertex;
+
+	/**
+	 * @brief Vertex with a point attached.
+	 */
+	class Simple_geometric_vertex : public Simple_vertex{
+		template<class ComplexGeometricTraits> friend class Skeleton_blocker_geometric_complex;
+	private:
+		Point point_;
+		Point& point(){	return point_; }
+		const Point& point() const {	return point_; }
+	};
+
+
+	class Simple_geometric_edge : public Skeleton_blocker_simple_traits::Graph_edge{
+		int index_;
+	public:
+		Simple_geometric_edge():Skeleton_blocker_simple_traits::Graph_edge(),index_(-1){}
+		/**
+		 * @brief Allows to modify the index of the edge.
+		 * The indices of the edge are used to store heap information
+		 * in the edge contraction algorithm.
+		 */
+		int& index(){return index_;}
+		int index() const {return index_;}
+	};
+
+
+	typedef Simple_geometric_vertex Graph_vertex;
+	typedef Skeleton_blocker_simple_traits::Graph_edge Graph_edge;
+};
+
+
+}
+
+}  // namespace GUDHI
+
+#endif /* GUDHI_SKELETON_BLOCKERS_SIMPLE_GEOMETRIC_TRAITS_H_ */
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h
new file mode 100644
index 00000000..3975bb46
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h
@@ -0,0 +1,142 @@
+/*
+ * Skeleton_blocker_simple_traits.h
+ *
+ *  Created on: Feb 11, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKERS_SIMPLE_TRAITS_H_
+#define GUDHI_SKELETON_BLOCKERS_SIMPLE_TRAITS_H_
+
+#include <string>
+#include <sstream>
+#include "Skeleton_blocker_simplex.h"
+
+namespace Gudhi{
+
+namespace skbl {
+
+/**
+ * @extends SkeletonBlockerDS
+ */
+struct Skeleton_blocker_simple_traits{
+	/**
+	 * @brief global and local handle similar to <a href="http://www.boost.org/doc/libs/1_38_0/libs/graph/doc/subgraph.html">boost subgraphs</a>.
+	 * In gross, vertices are stored in a vector.
+	 * For the root simplicial complex, the local and global descriptors are the same.
+	 * For a subcomplex L and one of its vertices 'v', the local descriptor of 'v' is its position in
+	 * the vertex vector of the subcomplex L whereas its global descriptor is the position of 'v'
+	 * in the vertex vector of the root simplicial complex.
+	 */
+	struct Root_vertex_handle{
+		typedef int boost_vertex_handle;
+		explicit Root_vertex_handle(boost_vertex_handle val = -1 ):vertex(val){}
+		boost_vertex_handle vertex;
+
+		bool operator!=( const Root_vertex_handle& other) const{
+			return ! (this->vertex == other.vertex);
+		}
+
+		bool operator==( const Root_vertex_handle& other) const{
+			return this->vertex == other.vertex;
+		}
+
+		bool operator<( const Root_vertex_handle& other) const{
+			return this->vertex < other.vertex;
+		}
+
+		friend std::ostream& operator << (std::ostream& o, const Root_vertex_handle & v)
+		{
+			o << v.vertex;
+			return o;
+		}
+	};
+
+	struct Vertex_handle{
+		typedef int boost_vertex_handle;
+		Vertex_handle(boost_vertex_handle val=-1):vertex(val){}
+
+		boost_vertex_handle vertex;
+
+		bool operator==( const Vertex_handle& other) const{
+			return this->vertex == other.vertex;
+		}
+
+		bool operator!=( const Vertex_handle& other) const{
+			return this->vertex != other.vertex;
+		}
+
+		bool operator<(const Vertex_handle& other) const{
+			return this->vertex < other.vertex;
+		}
+
+		friend std::ostream& operator << (std::ostream& o, const Vertex_handle & v)
+		{
+			o << v.vertex;
+			return o;
+		}
+	};
+
+
+
+	class Graph_vertex {
+		bool is_active_;
+		Root_vertex_handle id_;
+	public:
+		virtual ~Graph_vertex(){}
+
+		void activate(){is_active_=true;}
+		void deactivate(){is_active_=false;}
+		bool is_active() const{return is_active_;}
+		void set_id(Root_vertex_handle i){id_=i;}
+		Root_vertex_handle get_id() const{return id_;}
+
+		virtual std::string to_string() const {
+			std::ostringstream res;
+			res<< id_;
+			return res.str();
+		}
+
+		friend std::ostream& operator << (std::ostream& o, const Graph_vertex & v){
+			o << v.to_string();
+			return o;
+		}
+	};
+
+	class Graph_edge {
+		Root_vertex_handle a_;
+		Root_vertex_handle b_;
+		int index_;
+	public:
+
+		Graph_edge():a_(-1),b_(-1),index_(-1){}
+
+		int& index(){return index_;}
+		int index() const {return index_;}
+
+		void setId(Root_vertex_handle a,Root_vertex_handle b){
+			a_ = a;
+			b_ = b;
+		}
+
+		Root_vertex_handle first() const {
+			return a_;
+		}
+
+		Root_vertex_handle second() const {
+			return b_;
+		}
+
+		friend std::ostream& operator << (std::ostream& o, const Graph_edge & v){
+			o << "("<<v.a_<<","<<v.b_<<" - id = "<<v.index();
+			return o;
+		}
+	};
+
+};
+
+}
+
+}  // namespace GUDHI
+
+#endif /* GUDHI_SKELETON_BLOCKERS_SIMPLE_TRAITS_H_ */
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h
new file mode 100644
index 00000000..5593a5a8
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h
@@ -0,0 +1,375 @@
+#ifndef GUDHI_SKELETON_BLOCKER_SIMPLEX_H
+#define GUDHI_SKELETON_BLOCKER_SIMPLEX_H
+
+#include<cassert>
+#include<iostream>
+#include<set>
+#include<vector>
+#include <initializer_list>
+
+namespace Gudhi{
+
+namespace skbl {
+
+/**
+ *@brief Abstract simplex used in Skeleton blockers data-structure.
+ *
+ * An abstract simplex is represented as an ordered set of T elements,
+ * each element representing a vertex.
+ * 
+ * The element representing a vertex can be SkeletonBlockerDS::Vertex_handle or SkeletonBlockerDS::Root_vertex_handle.
+ *
+ *
+ */
+template <typename T>
+
+class Skeleton_blocker_simplex {
+
+private :
+	std::set<T> simplex_set;
+
+public:
+	typedef typename T::boost_vertex_handle boost_vertex_handle;
+
+	typedef T Vertex_handle;
+
+
+	typedef typename std::set<T>::const_iterator Simplex_vertex_const_iterator;
+	typedef typename std::set<T>::iterator Simplex_vertex_iterator;
+
+	/** @name Constructors / Destructors / Initialization
+	 */
+	//@{
+
+	
+//	Skeleton_blocker_simplex():simplex_set() {}
+
+	/**
+	 * Clear the simplex
+	 */
+	void clear() {
+		simplex_set.clear();
+	}
+
+
+	Skeleton_blocker_simplex(std::initializer_list<T>& list) {
+		for_each(list.begin(),list.end(),add_vertex);
+	}
+
+
+	template<typename... Args>
+	Skeleton_blocker_simplex(Args... args){
+		add_vertices(args...);
+	}
+
+
+	template<typename... Args>
+	void add_vertices(T v,Args... args){
+		add_vertex(v);
+		add_vertices(args...);
+	}
+
+	void add_vertices(T v){
+		add_vertex(v);
+	}
+
+	void add_vertices(){	
+	}
+
+	/**
+	 * Initialize a simplex with a string such as {0,1,2}
+	 */
+	Skeleton_blocker_simplex(std::string token){
+		clear();
+		if ((token[0] == '{')  && (token[token.size()-1] == '}' ) )
+		{
+			token.erase(0,1);
+			token.erase(token.size()-1,1);
+			while(token.size()!=0 ){
+				int coma_position=token.find_first_of(',');
+				//cout << "coma_position:"<<coma_position<<endl;
+				std::string n = token.substr(0,coma_position);
+				//cout << "token:"<<token<<endl;
+				token.erase(0,n.size()+1);
+				add_vertex((T)(atoi(n.c_str())));
+			}
+		}
+	}
+
+	//@}
+
+	/** @name Simplex manipulation
+	 */
+	//@{
+
+	/**
+	 * Add the vertex v to the simplex:
+	 * Note that adding two times the same vertex is
+	 * the same that adding it once.
+	 * \f$ (*this) \leftarrow (*this) \cup \{ v \} \f$
+	 */
+	void add_vertex(T v)
+	{
+		simplex_set.insert(v);
+	}
+
+	/**
+	 * Remove the vertex v from the simplex:
+	 * \f$ (*this) \leftarrow (*this) \setminus \{ v \} \f$
+	 */
+	void remove_vertex(T v)
+	{
+		simplex_set.erase(v);
+	}
+
+	/**
+	 * Intersects the simplex with the simplex a:
+	 * \f$ (*this) \leftarrow (*this) \cap a \f$
+	 */
+	void intersection(const Skeleton_blocker_simplex & a){
+		std::vector<T> v;
+		v.reserve(std::min(simplex_set.size(), a.simplex_set.size()));
+
+		set_intersection(simplex_set.begin(),simplex_set.end(),
+				a.simplex_set.begin(),a.simplex_set.end(),
+				std::back_inserter(v));
+		clear();
+		for (auto i:v)
+			simplex_set.insert(i);
+	}
+
+	/**
+	 * Substracts a from the simplex:
+	 * \f$ (*this) \leftarrow (*this) \setminus a \f$
+	 */
+	void difference(const Skeleton_blocker_simplex & a){
+		std::vector<T> v;
+		v.reserve(simplex_set.size());
+
+		set_difference(simplex_set.begin(),simplex_set.end(),
+				a.simplex_set.begin(),a.simplex_set.end(),
+				std::back_inserter(v));
+
+		clear();
+		for (auto i:v)
+			simplex_set.insert(i);
+	}
+
+	/**
+	 * Add vertices of a to the simplex:
+	 * \f$ (*this) \leftarrow (*this) \cup a \f$
+	 */
+	void union_vertices(const Skeleton_blocker_simplex & a){
+		std::vector<T> v;
+		v.reserve(simplex_set.size() + a.simplex_set.size());
+
+		set_union(simplex_set.begin(),simplex_set.end(),
+				a.simplex_set.begin(),a.simplex_set.end(),
+				std::back_inserter(v));
+		clear();
+		simplex_set.insert(v.begin(),v.end());
+	}
+
+	typename std::set<T>::const_iterator begin() const{
+		return simplex_set.cbegin();
+	}
+
+	typename std::set<T>::const_iterator end() const{
+		return simplex_set.cend();
+	}
+
+	typename std::set<T>::iterator begin(){
+		return simplex_set.begin();
+	}
+
+	typename std::set<T>::iterator end(){
+		return simplex_set.end();
+	}
+
+
+
+	//@}
+
+	/** @name Queries
+	 */
+	//@{
+
+	/**
+	 * Returns the dimension of the simplex.
+	 */
+	int dimension() const{
+		return (simplex_set.size() - 1);
+	}
+
+	bool empty() const{
+		return simplex_set.empty();
+	}
+
+	/**
+	 * Returns the first vertex of the (oriented) simplex.
+	 *
+	 * Be careful : assumes the simplex is non-empty.
+	 */
+	T first_vertex() const{
+		assert(!empty());
+		return *(begin());
+	}
+
+	/**
+	 * Returns the last vertex of the (oriented) simplex.
+	 *
+	 * Be careful : assumes the simplex is non-empty.
+	 */
+	T last_vertex() const
+	{
+		assert(!empty());
+		return *(simplex_set.rbegin());
+	}
+	/**
+	 * @return true iff the simplex contains the simplex a, i.e. iff \f$ a \subset (*this) \f$.
+	 */
+	bool contains(const Skeleton_blocker_simplex & a) const{
+		return includes(simplex_set.cbegin(),simplex_set.cend(),a.simplex_set.cbegin(),a.simplex_set.cend());
+	}
+
+	/**
+	 * @return true iff the simplex contains the difference \f$ a \setminus b \f$.
+	 */
+	bool contains_difference(const Skeleton_blocker_simplex& a, const Skeleton_blocker_simplex& b) const{
+		auto first1 = begin();
+		auto last1 = end();
+
+		auto first2 = a.begin();
+		auto last2 = a.end();
+
+		while (first2!=last2) {
+			// we ignore vertices of b
+			if(b.contains(*first2)){
+				++first2;
+			}
+			else{
+				if ( (first1==last1) || (*first2<*first1) ) return false;
+				if (!(*first1<*first2)) ++first2;
+				++first1;
+			}
+		}
+		return true;
+	}
+
+	/**
+	 * @return true iff the simplex contains the difference \f$ a \setminus \{ x \} \f$.
+	 */
+	bool contains_difference(const Skeleton_blocker_simplex& a, T x) const{
+		auto first1 = begin();
+		auto last1 = end();
+
+		auto first2 = a.begin();
+		auto last2 = a.end();
+
+		while (first2!=last2) {
+			// we ignore vertices of b
+			if(x == *first2){
+				++first2;
+			}
+			else{
+				if ( (first1==last1) || (*first2<*first1) ) return false;
+				if (!(*first1<*first2)) ++first2;
+				++first1;
+			}
+		}
+		return true;
+	}
+
+	/**
+	 * @return true iff the simplex contains the difference \f$ a \setminus \{ x \} \f$.
+	 */
+	bool contains_difference(const Skeleton_blocker_simplex& a, T x, T y) const{
+		auto first1 = begin();
+		auto last1 = end();
+
+		auto first2 = a.begin();
+		auto last2 = a.end();
+
+		while (first2!=last2) {
+			// we ignore vertices of b
+			if(x == *first2 || y == *first2){
+				++first2;
+			}
+			else{
+				if ( (first1==last1) || (*first2<*first1) ) return false;
+				if (!(*first1<*first2)) ++first2;
+				++first1;
+			}
+		}
+		return true;
+	}
+
+
+	/**
+	 * @return true iff the simplex contains the vertex v, i.e. iff \f$ v \in (*this) \f$.
+	 */
+	bool contains(T v) const{
+		return (simplex_set.find(v) != simplex_set.end());
+	}
+
+	/**
+	 * @return \f$ (*this) \cap a = \emptyset \f$.
+	 */
+	bool disjoint(const Skeleton_blocker_simplex& a) const{
+		std::vector<T> v;
+		v.reserve(std::min(simplex_set.size(), a.simplex_set.size()));
+
+		set_intersection(simplex_set.cbegin(),simplex_set.cend(),
+				a.simplex_set.cbegin(),a.simplex_set.cend(),
+				std::back_inserter(v));
+
+		return (v.size()==0);
+	}
+
+
+	bool operator==(const Skeleton_blocker_simplex& other) const{
+		return (this->simplex_set == other.simplex_set);
+	}
+
+	bool operator!=(const Skeleton_blocker_simplex& other) const{
+		return (this->simplex_set != other.simplex_set);
+	}
+
+	bool operator<(const Skeleton_blocker_simplex& other) const{
+		return (std::lexicographical_compare(this->simplex_set.begin(),this->simplex_set.end(),
+				other.begin(),other.end()));
+	}
+
+	//@}
+
+
+
+
+
+	/**
+	 * Display a simplex
+	 */
+	friend std::ostream& operator << (std::ostream& o, const Skeleton_blocker_simplex & sigma)
+	{
+		bool first = true;
+		o << "{";
+		for(auto i : sigma)
+		{
+			if(first) first = false ;
+			else o << ",";
+			o << i;
+		}
+		o << "}";
+		return o;
+	}
+
+
+};
+
+}
+
+}  // namespace GUDHI
+
+#endif
+
+
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h
new file mode 100644
index 00000000..5ca8e9af
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h
@@ -0,0 +1,295 @@
+#ifndef GUDHI_SKELETON_BLOCKER_SUB_COMPLEX_H
+#define GUDHI_SKELETON_BLOCKER_SUB_COMPLEX_H
+
+#include "gudhi/Skeleton_blocker_complex.h"
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h"
+#include "gudhi/Utils.h"
+
+namespace Gudhi{
+
+namespace skbl {
+
+/**
+ * @brief Simplicial subcomplex of a complex represented by a skeleton/blockers pair.
+ *
+ * @details Stores a subcomplex of a simplicial complex.
+ * To simplify explanations below, we will suppose that :
+ * - K is the root simplicial complex
+ * - L is a subcomplex of K.
+ *
+ * One vertex of K may exists in L but with a different address.
+ * To be able to locate the vertices in K from vertices of L, the class
+ * stores a map 'adresses' between vertices of K and vertices of L.
+ *
+ * Note that the type for handle of vertices of L is 'Vertex_handle' and
+ * the type for handle of vertices of K is 'Root_vertex_handle'.
+ *
+ * The template ComplexType is type of the root complex. It allows to know
+ * if the subcomplex is geometric or not.
+ * It has to be either 'Skeleton_blockers_complex' or 'Skeleton_blockers_geometric_complex'.
+ *
+ */
+template<typename ComplexType>
+class Skeleton_blocker_sub_complex : public ComplexType
+{
+
+protected:
+
+	template<class T> friend class Skeleton_blocker_link_complex;
+
+	typedef typename ComplexType::Graph Graph;
+	typedef typename ComplexType::Edge_handle Edge_handle;
+
+	typedef typename ComplexType::boost_vertex_handle boost_vertex_handle;
+
+
+public:
+	using ComplexType::add_vertex;
+	using ComplexType::add_edge;
+	using ComplexType::add_blocker;
+
+	typedef typename ComplexType::Vertex_handle Vertex_handle;
+	typedef typename ComplexType::Root_vertex_handle Root_vertex_handle;
+	typedef typename ComplexType::Simplex_handle Simplex_handle;
+	typedef typename ComplexType::Root_simplex_handle Root_simplex_handle;
+
+
+protected:
+
+
+	///**
+	//* @brief Returns true iff the simplex formed by all vertices contained in 'addresses_sigma_in_link'
+	//* but 'vertex_to_be_ignored' is in 'link'
+	//*/
+	/*
+			template<typename T> friend bool
+				proper_face_in_union(
+				Skeleton_blocker_sub_complex<T> & link,
+				std::vector<boost::optional<typename T::Vertex_handle> > & addresses_sigma_in_link,
+				int vertex_to_be_ignored);*/
+
+	/**
+	 * @brief Determines whether all proper faces of simplex 'sigma' belong to 'link1' \cup 'link2'
+	 * where 'link1' and 'link2' are subcomplexes of the same complex of type ComplexType
+	 */
+	//	template<typename T> friend	bool
+	//	proper_faces_in_union(Skeleton_blocker_simplex<typename T::Root_vertex_handle> & sigma, Skeleton_blocker_sub_complex<T> & link1, Skeleton_blocker_sub_complex<T> & link2){
+
+	//template<typename T> friend bool
+	//proper_faces_in_union(Skeleton_blocker_simplex<typename T::Root_vertex_handle> & sigma, Skeleton_blocker_sub_complex<T> & link1, Skeleton_blocker_sub_complex<T> & link2);
+
+	typedef std::map<Root_vertex_handle, Vertex_handle> IdAddressMap;
+	typedef typename IdAddressMap::value_type AddressPair;
+	typedef typename IdAddressMap::iterator IdAddressMapIterator;
+	typedef typename IdAddressMap::const_iterator IdAddressMapConstIterator;
+	std::map<Root_vertex_handle, Vertex_handle> adresses;
+
+
+public:
+	/**
+	 * Add a vertex 'global' of K to L. When added to L, this vertex will receive
+	 * another number, addresses(global), its local adress.
+	 * return the adress where the vertex lay on L.
+	 * The vertex corresponding to 'global' must not be already present
+	 * in the complex.
+	 */
+	Vertex_handle add_vertex(Root_vertex_handle global){
+		assert(!this->contains_vertex(global));
+		Vertex_handle address(boost::add_vertex(this->skeleton));
+		this->num_vertices_++;
+		(*this)[address].activate();
+		(*this)[address].set_id(global);
+		adresses.insert(AddressPair(global, address));
+		this->degree_.push_back(0);
+		return address;
+	}
+
+
+	/**
+	 * Add an edge (v1_root,v2_root) to the sub-complex.
+	 * It assumes that both vertices corresponding to v1_root and v2_root are present
+	 * in the sub-complex.
+	 */
+	void add_edge(Root_vertex_handle v1_root, Root_vertex_handle v2_root){
+		auto v1_sub(this->get_address(v1_root));
+		auto v2_sub(this->get_address(v2_root));
+		assert(v1_sub && v2_sub);
+		this->ComplexType::add_edge(*v1_sub, *v2_sub);
+	}
+
+	/**
+	 * Add a blocker to the sub-complex.
+	 * It assumes that all vertices of blocker_root are present
+	 * in the sub-complex.
+	 */
+	void add_blocker(const Root_simplex_handle& blocker_root){
+		auto blocker_sub = this->get_address(blocker_root);
+		assert(blocker_sub);
+		this->add_blocker(new Simplex_handle(*blocker_sub));
+	}
+
+
+
+public:
+
+	/**
+	 * Constructs the restricted complex of 'parent_complex' to
+	 * vertices of 'simplex'.
+	 */
+	friend void make_restricted_complex(const ComplexType & parent_complex, const Simplex_handle& simplex, Skeleton_blocker_sub_complex & result){
+		result.clear();
+		// add vertices to the sub complex
+		for (auto x : simplex){
+			assert(parent_complex.contains_vertex(x));
+			auto x_local = result.add_vertex(parent_complex[x].get_id());
+			result[x_local] = parent_complex[x];
+		}
+
+		// add edges to the sub complex
+		for (auto x : simplex){
+			// x_neigh is the neighbor of x intersected with vertices_simplex
+			Simplex_handle x_neigh;
+			parent_complex.add_neighbours(x, x_neigh, true);
+			x_neigh.intersection(simplex);
+			for (auto y : x_neigh){
+				result.add_edge(parent_complex[x].get_id(), parent_complex[y].get_id());
+			}
+		}
+
+		// add blockers to the sub complex
+		for (auto blocker : parent_complex.const_blocker_range()){
+			// check if it is the first time we encounter the blocker
+			if (simplex.contains(*blocker)){
+				Root_simplex_handle blocker_root(parent_complex.get_id(*(blocker)));
+				Simplex_handle blocker_restr(*result.get_simplex_address(blocker_root));
+				result.add_blocker(new Simplex_handle(blocker_restr));
+			}
+		}
+	}
+
+
+
+	void clear(){
+		adresses.clear();
+		ComplexType::clear();
+	}
+
+	/**
+	 * Compute the local vertex in L corresponding to the vertex global in K.
+	 * runs in O(log n) if n = num_vertices()
+	 */
+	boost::optional<Vertex_handle> get_address(Root_vertex_handle global) const{
+		boost::optional<Vertex_handle> res;
+		IdAddressMapConstIterator it = adresses.find(global);
+		if (it == adresses.end()) res.reset();
+		else  res = (*it).second;
+		return res;
+	}
+
+	//	/**
+	//	 * Allocates a simplex in L corresponding to the simplex s in K
+	//	 * with its local adresses and returns an AddressSimplex.
+	//	 */
+	//	boost::optional<Simplex_handle> get_address(const Root_simplex_handle & s) const;
+
+
+//private:
+	/**
+	 *  same as get_address except that it will return a simplex in any case.
+	 *  The vertices that were not found are not added.
+	 */
+	// @remark should be private but problem with VS
+	std::vector<boost::optional<Vertex_handle> > get_addresses(const Root_simplex_handle & s) const{
+		std::vector<boost::optional<Vertex_handle> > res;
+		for (auto i : s)
+		{
+			res.push_back(get_address(i));
+		}
+		return res;
+	}
+
+};
+
+
+/**
+ * @remark remarque perte de temps a creer un nouveau simplexe a chaque fois
+ * alors qu'on pourrait utiliser a la place de 'addresses_sigma_in_link'
+ * un simplex avec des valeurs spéciales ComplexDS::null_vertex par exemple
+ * pour indiquer qu'un vertex n'appartient pas au complex
+ */
+template<typename ComplexType>
+bool proper_face_in_union(
+		Skeleton_blocker_sub_complex<ComplexType> & link,
+		std::vector<boost::optional<typename ComplexType::Vertex_handle> > & addresses_sigma_in_link,
+		int vertex_to_be_ignored)
+{
+	// we test that all vertices of 'addresses_sigma_in_link' but 'vertex_to_be_ignored'
+	// are in link1 if it is the case we construct the corresponding simplex
+	bool vertices_sigma_are_in_link = true;
+	typename ComplexType::Simplex_handle sigma_in_link;
+	for (int i = 0; i < addresses_sigma_in_link.size(); ++i){
+		if (i != vertex_to_be_ignored){
+			if (!addresses_sigma_in_link[i]){
+				vertices_sigma_are_in_link = false;
+				break;
+			}
+			else sigma_in_link.add_vertex(*addresses_sigma_in_link[i]);
+		}
+	}
+	// If one of vertices of the simplex is not in the complex then it returns false
+	// Otherwise, it tests if the simplex is in the complex
+	return vertices_sigma_are_in_link && link.contains(sigma_in_link);
+}
+
+/*
+		template<typename ComplexType>
+		bool
+			proper_faces_in_union(Skeleton_blocker_simplex<typename ComplexType::Root_vertex_handle> & sigma, Skeleton_blocker_sub_complex<ComplexType> & link1, Skeleton_blocker_sub_complex<ComplexType> & link2)
+		{
+				typedef typename ComplexType::Vertex_handle  Vertex_handle;
+				std::vector<boost::optional<Vertex_handle> > addresses_sigma_in_link1 = link1.get_addresses(sigma);
+				std::vector<boost::optional<Vertex_handle> > addresses_sigma_in_link2 = link2.get_addresses(sigma);
+
+				for (int current_index = 0; current_index < addresses_sigma_in_link1.size(); ++current_index)
+				{
+
+					if (!proper_face_in_union(link1, addresses_sigma_in_link1, current_index)
+						&& !proper_face_in_union(link2, addresses_sigma_in_link2, current_index)){
+						return false;
+					}
+				}
+				return true;
+			}*/
+
+
+
+// Remark: this function should be friend in order to leave get_adresses private
+// however doing so seemes currently not possible due to a visual studio bug c2668 
+// "the compiler does not support partial ordering of template functions as specified in the C++ Standard"
+// http://www.serkey.com/error-c2668-ambiguous-call-to-overloaded-function-bb45ft.html
+template<typename ComplexType>
+bool
+proper_faces_in_union(Skeleton_blocker_simplex<typename ComplexType::Root_vertex_handle> & sigma, Skeleton_blocker_sub_complex<ComplexType> & link1, Skeleton_blocker_sub_complex<ComplexType> & link2)
+{
+	typedef typename ComplexType::Vertex_handle  Vertex_handle;
+	std::vector<boost::optional<Vertex_handle> > addresses_sigma_in_link1 = link1.get_addresses(sigma);
+	std::vector<boost::optional<Vertex_handle> > addresses_sigma_in_link2 = link2.get_addresses(sigma);
+
+	for (int current_index = 0; current_index < addresses_sigma_in_link1.size(); ++current_index)
+	{
+
+		if (!proper_face_in_union(link1, addresses_sigma_in_link1, current_index)
+				&& !proper_face_in_union(link2, addresses_sigma_in_link2, current_index)){
+			return false;
+		}
+	}
+	return true;
+}
+
+} // namespace skbl
+
+}  // namespace GUDHI
+
+
+#endif
+
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Top_faces.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Top_faces.h
new file mode 100644
index 00000000..0d870f1a
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Top_faces.h
@@ -0,0 +1,53 @@
+/*
+ * Top_faces.h
+ *
+ *  Created on: Oct 23, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef TOP_FACES_H_
+#define TOP_FACES_H_
+
+#include <list>
+#include <vector>
+#include <set>
+
+template<typename SimplexHandle>
+std::list<SimplexHandle> subfaces(SimplexHandle top_face){
+	std::list<SimplexHandle> res;
+	if(top_face.dimension()==-1) return res;
+	if(top_face.dimension()==0) {
+		res.push_back(top_face);
+		return res;
+	}
+	else{
+		auto first_vertex = top_face.first_vertex();
+		top_face.remove_vertex(first_vertex);
+		res = subfaces(top_face);
+		std::list<SimplexHandle> copy = res;
+		for(auto& simplex : copy){
+			simplex.add_vertex(first_vertex);
+		}
+		res.push_back(SimplexHandle(first_vertex));
+		res.splice(res.end(),copy);
+		return res;
+	}
+}
+
+/**
+ * add all faces of top_face in simplices_per_dimension
+ */
+template<typename SimplexHandle>
+void register_faces(
+		std::vector< std::set<SimplexHandle> >& simplices_per_dimension,
+		const SimplexHandle& top_face){
+	std::list<SimplexHandle> subfaces_list = subfaces(top_face);
+	for(auto& simplex : subfaces_list ){
+		simplices_per_dimension[simplex.dimension()].insert(simplex);
+	}
+}
+
+
+
+
+#endif /* TOP_FACES_H_ */
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h
new file mode 100644
index 00000000..21091d10
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h
@@ -0,0 +1,120 @@
+/*
+ * Skeleton_blockers_blockers_iterators.h
+ *
+ *  Created on: Aug 25, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKERS_BLOCKERS_ITERATORS_H_
+#define GUDHI_SKELETON_BLOCKERS_BLOCKERS_ITERATORS_H_
+
+#include "boost/iterator/iterator_facade.hpp"
+
+namespace Gudhi{
+
+namespace skbl {
+
+/**
+ * @brief Iterator through the blockers of a vertex.
+  */
+// ReturnType = const Simplex_handle* or Simplex_handle*
+// MapIteratorType = BlockerMapConstIterator or BlockerMapIterator
+template<typename MapIteratorType, typename ReturnType>
+class Blocker_iterator_internal : public boost::iterator_facade<
+  Blocker_iterator_internal<MapIteratorType,ReturnType>,
+  ReturnType,
+  boost::forward_traversal_tag,
+  ReturnType
+  >{
+private:
+	MapIteratorType current_position;
+	MapIteratorType end_of_map;
+public:
+
+	Blocker_iterator_internal():current_position(){}
+
+	Blocker_iterator_internal(MapIteratorType position,MapIteratorType end_of_map_ ):
+		current_position(position), end_of_map(end_of_map_)
+	{	}
+
+	bool equal(const Blocker_iterator_internal& other) const{
+		return current_position == other.current_position;
+	}
+
+	void increment(){
+		goto_next_blocker();
+	}
+
+	ReturnType dereference() const	{
+		return(current_position->second);
+	}
+
+private:
+	/**
+	 * Let the current pair be (v,sigma) where v is a vertex and sigma is a blocker.
+	 * If v is not the first vertex of sigma then we already have seen sigma as a blocker
+	 * and we look for the next one.
+	 */
+	void goto_next_blocker(){
+		do {
+			++current_position;
+		} while (!(current_position == end_of_map) && !first_time_blocker_is_seen());
+	}
+
+	bool first_time_blocker_is_seen() const{
+		return current_position->first  == current_position->second->first_vertex();
+	}
+};
+
+
+
+/**
+ * @brief Iterator through the blockers of a vertex
+ */
+// ReturnType = const Simplex_handle* or Simplex_handle*
+// MapIteratorType = BlockerMapConstIterator or BlockerMapIterator
+template<typename MapIteratorType, typename ReturnType>
+class Blocker_iterator_around_vertex_internal : public boost::iterator_facade<
+	Blocker_iterator_around_vertex_internal<MapIteratorType,ReturnType>,
+	ReturnType,
+	boost::forward_traversal_tag,
+	ReturnType
+>{
+private:
+	MapIteratorType current_position_;
+public:
+
+	Blocker_iterator_around_vertex_internal():current_position_(){}
+
+	Blocker_iterator_around_vertex_internal(MapIteratorType position):
+		current_position_(position)
+	{}
+
+	Blocker_iterator_around_vertex_internal& operator=(Blocker_iterator_around_vertex_internal other){
+		this->current_position_ = other.current_position_;
+		return *this;
+	}
+
+	bool equal(const Blocker_iterator_around_vertex_internal& other) const{
+		return current_position_ == other.current_position_;
+	}
+
+	void increment(){
+		current_position_++;
+	}
+
+	ReturnType dereference() const{
+		return(current_position_->second);
+	}
+
+
+	MapIteratorType current_position(){
+		return this->current_position_;
+	}
+};
+
+}
+
+} // namespace GUDHI
+
+#endif /* GUDHI_SKELETON_BLOCKERS_BLOCKERS_ITERATORS_H_ */
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h
new file mode 100644
index 00000000..82025a1c
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h
@@ -0,0 +1,153 @@
+/*
+ * Skeleton_blockers_iterators_out.h
+ *
+ *  Created on: Aug 25, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKERS_ITERATORS_EDGES_H_
+#define GUDHI_SKELETON_BLOCKERS_ITERATORS_EDGES_H_
+
+#include "boost/iterator/iterator_facade.hpp"
+
+
+namespace Gudhi{
+
+namespace skbl {
+
+template<typename SkeletonBlockerComplex>
+class Complex_edge_around_vertex_iterator :
+ public boost::iterator_facade < Complex_edge_around_vertex_iterator<SkeletonBlockerComplex>
+		, typename SkeletonBlockerComplex::Edge_handle const
+		, boost::forward_traversal_tag
+		, typename SkeletonBlockerComplex::Edge_handle const
+		>
+{
+	friend class boost::iterator_core_access;
+
+	typedef SkeletonBlockerComplex Complex;
+	typedef typename Complex::boost_adjacency_iterator boost_adjacency_iterator;
+	typedef typename Complex::Vertex_handle Vertex_handle;
+	typedef typename Complex::Edge_handle Edge_handle;
+
+private:
+
+	const Complex* complex;
+	Vertex_handle v;
+
+	boost_adjacency_iterator current_;
+	boost_adjacency_iterator end_;
+
+public:
+
+	Complex_edge_around_vertex_iterator():complex(NULL){
+	}
+
+	Complex_edge_around_vertex_iterator(const Complex* complex_,Vertex_handle v_):
+		complex(complex_),
+		v(v_)
+	{
+		tie(current_,end_) = adjacent_vertices(v.vertex, complex->skeleton);
+	}
+
+	/**
+	 * returns an iterator to the end
+	 */
+	Complex_edge_around_vertex_iterator(const Complex* complex_,Vertex_handle v_,int end):
+		complex(complex_),
+		v(v_)
+	{
+		tie(current_,end_) = adjacent_vertices(v.vertex, complex->skeleton);
+		set_end();
+	}
+
+	bool equal(const Complex_edge_around_vertex_iterator& other) const{
+		return (complex== other.complex) && (v == other.v) && (current_ == other.current_) && (end_ == other.end_);
+	}
+
+	void increment(){
+		if(current_ != end_)
+			++current_;
+	}
+
+	Edge_handle dereference() const{
+		return *(*complex)[std::make_pair(v,*current_)];
+	}
+
+private:
+	//remove this ugly hack
+	void set_end(){
+		current_ = end_;
+	}
+};
+
+
+
+/**
+ *@brief Iterator on the edges of a simplicial complex.
+ *
+ */
+template<typename SkeletonBlockerComplex>
+class Complex_edge_iterator :
+public boost::iterator_facade < Complex_edge_iterator<SkeletonBlockerComplex>
+, typename SkeletonBlockerComplex::Edge_handle const
+, boost::forward_traversal_tag
+, typename SkeletonBlockerComplex::Edge_handle const
+>
+
+{
+	friend class boost::iterator_core_access;
+public:
+	typedef SkeletonBlockerComplex Complex;
+	typedef typename Complex::boost_edge_iterator boost_edge_iterator;
+	typedef typename Complex::Edge_handle Edge_handle;
+
+
+	const Complex* complex;
+	std::pair<boost_edge_iterator,boost_edge_iterator> edge_iterator ;
+
+	Complex_edge_iterator():complex(NULL){
+	}
+
+	Complex_edge_iterator(const SkeletonBlockerComplex* complex_):
+		complex(complex_),
+		edge_iterator(boost::edges(complex_->skeleton))
+	{
+	}
+
+	/**
+	 * return an iterator to the end
+	 */
+	Complex_edge_iterator(const SkeletonBlockerComplex* complex_,int end):
+		complex(complex_),
+		edge_iterator(boost::edges(complex_->skeleton))
+	{
+		edge_iterator.first = edge_iterator.second;
+	}
+
+
+	bool equal(const Complex_edge_iterator& other) const{
+		return (complex == other.complex) && (edge_iterator == other.edge_iterator);
+	}
+
+	void increment(){
+		if(edge_iterator.first != edge_iterator.second){
+			++(edge_iterator.first);
+		}
+	}
+
+	Edge_handle dereference() const{
+		return(*(edge_iterator.first));
+	}
+};
+
+
+
+}
+
+} // namespace GUDHI
+
+
+#endif /* GUDHI_SKELETON_BLOCKERS_ITERATORS_EDGES_H_ */
+
+
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h
new file mode 100644
index 00000000..cf827705
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h
@@ -0,0 +1,21 @@
+/*
+ * Skeleton_blockers_iterators.h
+ *
+ *  Created on: Aug 25, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKERS_ITERATORS_H_
+#define GUDHI_SKELETON_BLOCKERS_ITERATORS_H_
+
+
+#include "Skeleton_blockers_vertices_iterators.h"
+#include "Skeleton_blockers_edges_iterators.h"
+#include "Skeleton_blockers_blockers_iterators.h"
+#include "Skeleton_blockers_triangles_iterators.h"
+#include "Skeleton_blockers_simplices_iterators.h"
+
+
+
+
+#endif /* GUDHI_SKELETON_BLOCKERS_ITERATORS_H_ */
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h
new file mode 100644
index 00000000..5ac0f21e
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h
@@ -0,0 +1,412 @@
+/*
+ * Skeleton_blockers_simplices_iterators.h
+ *
+ *  Created on: Sep 26, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_KELETON_BLOCKERS_SIMPLICES_ITERATORS_H_
+#define GUDHI_SKELETON_BLOCKERS_SIMPLICES_ITERATORS_H_
+
+#include <memory>
+#include <list>
+#include <iostream>
+#include "gudhi/Utils.h"
+#include "boost/iterator/iterator_facade.hpp"
+
+
+#include "gudhi/Skeleton_blocker_link_complex.h"
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h"
+
+
+
+
+namespace Gudhi {
+
+
+namespace skbl {
+
+
+/**
+ * Link may be Skeleton_blocker_link_complex<SkeletonBlockerComplex> to iterate over all
+ * simplices around a vertex OR
+ * Skeleton_blocker_superior_link_complex<SkeletonBlockerComplex> to iterate over all
+ * superior vertices around a vertex.
+ * The iteration is done by computing a trie with the link and doing a breadth-first traversal
+ * of the trie.
+ */
+template<typename SkeletonBlockerComplex,typename Link>
+class Simplex_around_vertex_iterator :
+		public boost::iterator_facade < Simplex_around_vertex_iterator<SkeletonBlockerComplex,Link>
+, typename SkeletonBlockerComplex::Simplex_handle
+, boost::forward_traversal_tag
+, typename SkeletonBlockerComplex::Simplex_handle
+>
+{
+	friend class boost::iterator_core_access;
+	typedef SkeletonBlockerComplex Complex;
+	typedef typename Complex::Vertex_handle Vertex_handle;
+	typedef typename Complex::Edge_handle Edge_handle;
+	typedef typename Complex::Simplex_handle Simplex_handle;
+
+
+	typedef typename Link::Vertex_handle Link_vertex_handle;
+	// Link_vertex_handle == Complex_Vertex_handle but this renaming helps avoiding confusion
+
+
+private:
+
+	struct Trie{
+		Vertex_handle v;
+		std::vector<std::shared_ptr<Trie> > childs;
+	private:
+		const Trie* parent_;
+	public:
+
+
+		//std::list<std::unique_ptr<Trie> > childs; -> use of deleted function
+
+		Trie():parent_(0){}
+		Trie(Vertex_handle v_):v(v_),parent_(0){
+		}
+
+		Trie(Vertex_handle v_,Trie* parent):v(v_),parent_(parent){
+		}
+
+
+		bool operator==(const Trie& other) const{
+			return (v == other.v) ;
+		}
+
+		void add_child(Trie* child){
+			if(child){
+				std::shared_ptr<Trie> ptr_to_add(child);
+				childs.push_back(ptr_to_add);
+				child->parent_ = this;
+			}
+		}
+
+
+		friend std::ostream& operator<<(std::ostream& stream, const Trie& trie){
+			stream<< "T( "<< trie.v<< " ";
+			for(auto t : trie.childs)
+				stream << *t ;
+			stream<<")";
+			return stream;
+		}
+
+		// goes to the root in the trie to consitute simplex
+		void add_vertices_up_to_the_root(Simplex_handle& res) const{
+			res.add_vertex(v);
+			if(parent_)
+				parent_->add_vertices_up_to_the_root(res);
+		}
+
+		Simplex_handle simplex() const{
+			Simplex_handle res;
+			add_vertices_up_to_the_root(res);
+			return res;
+		}
+
+		bool is_leaf() const{
+			return childs.empty();
+		}
+
+		bool is_root() const{
+			return parent_==0;
+		}
+
+		const Trie* parent() {
+			return parent_;
+		}
+
+		void remove_leaf() {
+			assert(is_leaf);
+			if(!is_root())
+				parent_->childs.erase(this);
+		}
+
+	private:
+
+
+	public:
+
+		Trie* go_bottom_left(){
+			if(is_leaf())
+				return this;
+			else
+				return (*childs.begin())->go_bottom_left();
+		}
+
+	};
+
+private:
+	const Complex* complex;
+	Vertex_handle v;
+	std::shared_ptr<Link> link_v;
+	std::shared_ptr<Trie> trie;
+	std::list<Trie*> nodes_to_be_seen; // todo regrouper
+
+public:
+	Simplex_around_vertex_iterator():complex(0){
+	}
+
+	Simplex_around_vertex_iterator(const Complex* complex_,Vertex_handle v_):
+		complex(complex_),
+		v(v_),
+		link_v(new Link(*complex_,v_)),
+		trie(new Trie(v_)){
+		compute_trie_and_nodes_to_be_seen();
+	}
+
+	// todo avoid useless copy
+	// todo currently just work if copy begin iterator
+	Simplex_around_vertex_iterator(const Simplex_around_vertex_iterator& other):
+		complex(other.complex),
+		v(other.v),
+		link_v(other.link_v),
+		trie(other.trie),
+		nodes_to_be_seen(other.nodes_to_be_seen)
+	{
+		if(!other.is_end()){
+		}
+	}
+
+	/**
+	 * returns an iterator to the end
+	 */
+	Simplex_around_vertex_iterator(const Complex* complex_,Vertex_handle v_,bool end):
+		complex(complex_),
+		v(v_){
+		set_end();
+	}
+
+private:
+
+
+	void compute_trie_and_nodes_to_be_seen(){
+		// once we go through every simplices passing through v0
+		// we remove v0. That way, it prevents from passing a lot of times
+		// though edges leaving v0.
+		// another solution would have been to provides an adjacency iterator
+		// to superior vertices that avoids lower ones.
+		while(!link_v->empty()){
+			auto v0 = *(link_v->vertex_range().begin());
+			trie->add_child(build_trie(v0,trie.get()));
+			link_v->remove_vertex(v0);
+		}
+		nodes_to_be_seen.push_back(trie.get());
+	}
+
+	Trie* build_trie(Link_vertex_handle link_vh,Trie* parent){
+		Trie* res = new Trie(parent_vertex(link_vh),parent);
+		for(Link_vertex_handle nv : link_v->vertex_range(link_vh)) {
+			if(link_vh < nv){
+				Simplex_handle simplex_node_plus_nv(res->simplex());
+				simplex_node_plus_nv.add_vertex(parent_vertex(nv));
+				if(complex->contains(simplex_node_plus_nv)){
+					res->add_child(build_trie(nv,res));
+				}
+			}
+		}
+		return res;
+	}
+
+	bool is_node_in_complex(Trie* trie){
+		return true;
+	}
+
+	Vertex_handle parent_vertex(Link_vertex_handle link_vh) const{
+		return complex->convert_handle_from_another_complex(*link_v,link_vh);
+	}
+
+
+
+public:
+
+	friend std::ostream& operator<<(std::ostream& stream, const Simplex_around_vertex_iterator& savi){
+		stream<< savi.trie<< std::endl; ;
+		stream << "("<<savi.nodes_to_be_seen.size()<<") nodes to see\n";
+		return stream;
+	}
+
+	bool equal(const Simplex_around_vertex_iterator& other) const{
+		bool same_complex = (complex == other.complex);
+		if(!same_complex)
+			return false;
+
+		if(!(v == other.v))
+			return false;
+
+		bool both_empty = nodes_to_be_seen.empty() && other.nodes_to_be_seen.empty();
+		if(both_empty)
+			return true;
+
+		bool both_non_empty = !nodes_to_be_seen.empty() && !other.nodes_to_be_seen.empty();
+
+		if(!both_non_empty) return false; //one is empty the other is not
+
+		bool same_node = (**(nodes_to_be_seen.begin()) == **(other.nodes_to_be_seen.begin()));
+		return same_node;
+	}
+
+	void increment(){
+		assert(!is_end());
+		Trie* first_node = nodes_to_be_seen.front();
+
+		nodes_to_be_seen.pop_front();
+
+		for(auto childs : first_node->childs){
+			nodes_to_be_seen.push_back(childs.get());
+		}
+
+	}
+
+	Simplex_handle dereference() const{
+		assert(!nodes_to_be_seen.empty());
+		Trie* first_node = nodes_to_be_seen.front();
+		return first_node->simplex();
+	}
+
+private:
+	void set_end(){
+		nodes_to_be_seen.clear();
+	}
+
+	bool is_end() const{
+		return nodes_to_be_seen.empty();
+	}
+};
+
+
+
+template<typename SkeletonBlockerComplex>
+class Simplex_iterator :
+		public boost::iterator_facade < Simplex_iterator<SkeletonBlockerComplex>
+, typename SkeletonBlockerComplex::Simplex_handle
+, boost::forward_traversal_tag
+, typename SkeletonBlockerComplex::Simplex_handle
+>
+{
+	typedef Skeleton_blocker_link_superior<SkeletonBlockerComplex> Link;
+
+	friend class boost::iterator_core_access;
+
+	template<class SkBlDS> friend class Skeleton_blocker_complex;
+
+
+	typedef SkeletonBlockerComplex Complex;
+	typedef typename Complex::Vertex_handle Vertex_handle;
+	typedef typename Complex::Edge_handle Edge_handle;
+	typedef typename Complex::Simplex_handle Simplex_handle;
+
+	typedef typename Complex::CVI CVI;
+
+
+	typedef typename Link::Vertex_handle Link_vertex_handle;
+
+private:
+
+	const Complex* complex_;
+	CVI current_vertex_;
+
+	typedef Simplex_around_vertex_iterator<SkeletonBlockerComplex,Link> SAVI;
+	SAVI current_simplex_around_current_vertex_;
+	SAVI simplices_around_current_vertex_end_;
+
+
+public:
+	Simplex_iterator():complex_(0){}
+
+	// should not be called if the complex is empty
+	Simplex_iterator(const Complex* complex):
+		complex_(complex),
+		current_vertex_(complex->vertex_range().begin()),
+		current_simplex_around_current_vertex_(complex,*current_vertex_),
+		simplices_around_current_vertex_end_(complex,*current_vertex_,true)
+	{
+		assert(!complex->empty());
+	}
+
+private:
+	// todo return to private
+	Simplex_iterator(const Complex* complex,bool end):
+		complex_(complex)
+	{
+		set_end();
+	}
+
+public:
+
+	Simplex_iterator(const Simplex_iterator& other)
+	:
+		complex_(other.complex_),
+		current_vertex_(other.current_vertex_),
+		current_simplex_around_current_vertex_(other.current_simplex_around_current_vertex_),
+		simplices_around_current_vertex_end_(other.simplices_around_current_vertex_end_)
+	{
+	}
+
+	friend Simplex_iterator make_begin_iterator(const Complex* complex){
+		if(complex->empty())
+			return make_end_simplex_iterator(complex);
+		else
+			return Simplex_iterator(complex);
+	}
+
+	friend Simplex_iterator make_end_simplex_iterator(const Complex* complex){
+		return Simplex_iterator(complex,true);
+	}
+
+	bool equal(const Simplex_iterator& other) const{
+		if(complex_!=other.complex_) return false;
+		if(current_vertex_!=other.current_vertex_) return false;
+		if(is_end() && other.is_end()) return true;
+		if(current_simplex_around_current_vertex_ != other.current_simplex_around_current_vertex_)
+			return false;
+		return true;
+	}
+
+	void increment(){
+		if(current_simplex_around_current_vertex_!= simplices_around_current_vertex_end_){
+			current_simplex_around_current_vertex_.increment();
+			if( current_simplex_around_current_vertex_== simplices_around_current_vertex_end_)
+				goto_next_vertex();
+		}
+		else{
+			goto_next_vertex();
+		}
+	}
+
+	void goto_next_vertex(){
+		current_vertex_.increment();
+		if(!is_end()){
+			current_simplex_around_current_vertex_= SAVI(complex_,*current_vertex_);
+			simplices_around_current_vertex_end_ = SAVI(complex_,*current_vertex_,true);
+		}
+	}
+
+	Simplex_handle dereference() const{
+		return current_simplex_around_current_vertex_.dereference();
+	}
+
+private:
+	void set_end(){
+		current_vertex_ = complex_->vertex_range().end();
+	}
+
+	bool is_end() const{
+		return (current_vertex_ == complex_->vertex_range().end());
+	}
+};
+
+
+}
+
+} // namespace GUDHI
+
+
+
+
+
+#endif /* SKELETON_BLOCKERS_SIMPLICES_ITERATORS_H_ */
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h
new file mode 100644
index 00000000..4b66ced5
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h
@@ -0,0 +1,224 @@
+/*
+ * Skeleton_blockers_triangles_iterators.h
+ *
+ *  Created on: Aug 25, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKERS_TRIANGLES_ITERATORS_H_
+#define GUDHI_SKELETON_BLOCKERS_TRIANGLES_ITERATORS_H_
+
+#include "boost/iterator/iterator_facade.hpp"
+#include <memory>
+
+namespace Gudhi{
+
+namespace skbl {
+
+//////////////////////////////////////////////////////////////////////
+/**
+ * \brief Iterator over the triangles that are
+ * adjacent to a vertex of the simplicial complex.
+ * \remark Will be removed soon -> dont look
+ */
+template<typename Complex,typename LinkType>
+class Triangle_around_vertex_iterator : public boost::iterator_facade
+< Triangle_around_vertex_iterator <Complex,LinkType>
+, typename Complex::Simplex_handle const
+, boost::forward_traversal_tag
+, typename Complex::Simplex_handle const
+>
+{
+	friend class boost::iterator_core_access;
+	template<typename T> friend class Triangle_iterator;
+private:
+	typedef typename LinkType::Vertex_handle Vertex_handle;
+	typedef typename LinkType::Root_vertex_handle Root_vertex_handle;
+	typedef typename LinkType::Simplex_handle Simplex_handle;
+	typedef Complex_edge_iterator<Complex> Complex_edge_iterator_;
+
+	const Complex* complex_;
+	Vertex_handle v_;
+	std::shared_ptr<LinkType> link_;
+	Complex_edge_iterator_ current_edge_;
+	bool is_end_;
+public:
+	Triangle_around_vertex_iterator(const Complex* complex,Vertex_handle v):
+		complex_(complex),v_(v),link_(new LinkType(*complex,v_)),
+		current_edge_(link_->edge_range().begin()),
+		is_end_(current_edge_ == link_->edge_range().end()){
+	}
+
+	/**
+	 * @brief ugly hack to get an iterator to the end
+	 */
+	Triangle_around_vertex_iterator(const Complex* complex,Vertex_handle v,bool is_end):
+		complex_(complex),v_(v),link_(0),is_end_(true){
+	}
+
+	/**
+	 * @brief ugly hack to get an iterator to the end
+	 */
+	Triangle_around_vertex_iterator():
+		complex_(0),v_(-1),link_(0),is_end_(true){
+	}
+
+
+	Triangle_around_vertex_iterator(const Triangle_around_vertex_iterator& other){
+		v_ = other.v_;
+		complex_ = other.complex_;
+		is_end_ = other.is_end_;
+
+		if(!is_end_){
+			link_ = other.link_;
+			current_edge_= other.current_edge_;
+		}
+	}
+
+	bool equal(const Triangle_around_vertex_iterator& other) const{
+		return (complex_==other.complex_) && ((finished() &&other.finished()) || current_edge_ == other.current_edge_);
+	}
+
+	Simplex_handle dereference() const{
+		Root_vertex_handle v1 = (*link_)[*current_edge_].first();
+		Root_vertex_handle v2 = (*link_)[*current_edge_].second();
+		return Simplex_handle(v_,*(complex_->get_address(v1)),*(complex_->get_address(v2)));
+	}
+
+	void increment(){
+		++current_edge_;
+	}
+
+private:
+	bool finished() const{
+		return is_end_ || (current_edge_ == link_->edge_range().end());
+	}
+
+};
+
+
+
+/**
+ * \brief Iterator over the triangles of the
+ * simplicial complex.
+ * \remark Will be removed soon -> dont look
+ *
+ */
+template<typename SkeletonBlockerComplex>
+class Triangle_iterator :
+		public boost::iterator_facade<
+		Triangle_iterator <SkeletonBlockerComplex>,
+		typename SkeletonBlockerComplex::Simplex_handle const
+		, boost::forward_traversal_tag
+		, typename SkeletonBlockerComplex::Simplex_handle const
+		>
+{
+	friend class boost::iterator_core_access;
+private:
+	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::Superior_triangle_around_vertex_iterator STAVI;
+
+	const SkeletonBlockerComplex* complex_;
+	Complex_vertex_iterator<SkeletonBlockerComplex> current_vertex_;
+	STAVI current_triangle_;
+	bool is_end_;
+public:
+
+	/*
+	 * @remark  assume that the complex is non-empty
+	 */
+	Triangle_iterator(const SkeletonBlockerComplex* complex):
+		complex_(complex),
+		current_vertex_(complex->vertex_range().begin()),
+		current_triangle_(complex,*current_vertex_), // xxx this line is problematic is the complex is empty
+		is_end_(false){
+
+		assert(!complex->empty());
+		gotoFirstTriangle();
+	}
+
+private:
+	//goto to the first triangle or to the end if none
+	void gotoFirstTriangle(){
+		if(!is_finished() && current_triangle_.finished()){
+			goto_next_vertex();
+		}
+	}
+public:
+
+	/**
+	 * @brief ugly hack to get an iterator to the end
+	 * @remark  assume that the complex is non-empty
+	 */
+	Triangle_iterator(const SkeletonBlockerComplex* complex,bool is_end):
+		complex_(complex),
+		current_vertex_(complex->vertex_range().end()),
+		current_triangle_(), // xxx this line is problematic is the complex is empty
+		is_end_(true){
+	}
+
+
+	Triangle_iterator& operator=(const Triangle_iterator & other){
+		complex_ = other.complex_;
+		Complex_vertex_iterator<SkeletonBlockerComplex> current_vertex_;
+		STAVI current_triangle_;
+		return *this;
+	}
+
+
+	bool equal(const Triangle_iterator& other) const{
+		bool both_are_finished = is_finished() && other.is_finished();
+		bool both_arent_finished = !is_finished() && !other.is_finished();
+		// if the two iterators are not finished, they must have the same state
+		return (complex_==other.complex_) &&
+				(both_are_finished	||
+						( (both_arent_finished) && current_vertex_ == other.current_vertex_ && current_triangle_ == other.current_triangle_));
+
+	}
+
+	Simplex_handle dereference() const{
+		return *current_triangle_;
+	}
+
+private:
+
+	// goto the next vertex that has a triangle pending or the
+	// end vertex iterator if none exists
+	void goto_next_vertex(){
+		assert(current_triangle_.finished()); //we mush have consume all triangles passing through the vertex
+		assert(!is_finished()); // we must not be done
+
+		++current_vertex_;
+
+		if(!is_finished()){
+			current_triangle_ = STAVI(complex_, *current_vertex_);
+			if(current_triangle_.finished())
+				goto_next_vertex();
+		}
+	}
+public:
+	void increment(){
+		if(!current_triangle_.finished()){
+			++current_triangle_; // problem here
+			if(current_triangle_.finished())
+				goto_next_vertex();
+		}
+		else{
+			assert(!is_finished());
+			goto_next_vertex();
+		}
+	}
+
+private:
+	bool is_finished() const{
+		return is_end_ || current_vertex_ == complex_->vertex_range().end();
+	}
+};
+
+}
+
+} // namespace GUDHI
+
+#endif /* GUDHI_SKELETON_BLOCKERS_TRIANGLES_ITERATORS_H_ */
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h
new file mode 100644
index 00000000..875e915a
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h
@@ -0,0 +1,168 @@
+/*
+ * Skeleton_blockers_iterators_out.h
+ *
+ *  Created on: Aug 25, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKERS_VERTICES_ITERATORS_H_
+#define GUDHI_SKELETON_BLOCKERS_VERTICES_ITERATORS_H_
+
+#include "boost/iterator/iterator_facade.hpp"
+
+
+namespace Gudhi{
+
+namespace skbl {
+
+/**
+ *@brief Iterator on the vertices of a simplicial complex
+ *@remark The increment operator go to the next active vertex.
+ *@remark Incrementation increases Vertex_handle.
+ */
+template<typename SkeletonBlockerComplex>
+class Complex_vertex_iterator : public boost::iterator_facade
+< Complex_vertex_iterator <SkeletonBlockerComplex>
+  , typename SkeletonBlockerComplex::Vertex_handle const
+  , boost::forward_traversal_tag
+  , typename SkeletonBlockerComplex::Vertex_handle const
+  >
+{
+	friend class boost::iterator_core_access;
+
+	typedef typename SkeletonBlockerComplex::boost_vertex_iterator boost_vertex_iterator;
+	typedef typename SkeletonBlockerComplex::Vertex_handle Vertex_handle;
+private:
+	const SkeletonBlockerComplex* complex;
+	std::pair<boost_vertex_iterator,boost_vertex_iterator> vertexIterator;
+
+
+public:
+	Complex_vertex_iterator():complex(NULL){
+	}
+
+	Complex_vertex_iterator(const SkeletonBlockerComplex* complex_):
+		complex(complex_),
+		vertexIterator(vertices(complex_->skeleton)){
+		if(!finished() && !is_active()) {
+			goto_next_valid();
+		}
+	}
+
+	/**
+	 * return an iterator to the end.
+	 */
+	Complex_vertex_iterator(const SkeletonBlockerComplex* complex_,int end):
+		complex(complex_),vertexIterator(vertices(complex_->skeleton)){
+		vertexIterator.first = vertexIterator.second ;
+	}
+
+public:
+	void increment () {goto_next_valid();}
+	Vertex_handle  dereference() const	{
+		return(Vertex_handle(*(vertexIterator.first)));
+	}
+
+	bool equal(const Complex_vertex_iterator& other) const{
+		return vertexIterator == other.vertexIterator &&  complex == other.complex;
+	}
+
+	bool operator<(const Complex_vertex_iterator& other) const{
+		return dereference()<other.dereference();
+	}
+
+private:
+	bool finished() const{
+		return vertexIterator.first == vertexIterator.second;
+	}
+
+	void goto_next_valid(){
+		++vertexIterator.first;
+		if(!finished() && !is_active()){
+			goto_next_valid();
+		}
+	}
+
+	bool is_active() const{
+		return ((*complex)[Vertex_handle(*vertexIterator.first)]).is_active();
+	}
+
+};
+
+
+
+
+template<typename SkeletonBlockerComplex>
+class Complex_neighbors_vertices_iterator
+: public boost::iterator_facade < Complex_neighbors_vertices_iterator<SkeletonBlockerComplex>
+  , typename SkeletonBlockerComplex::Vertex_handle const
+  , boost::forward_traversal_tag
+  , typename SkeletonBlockerComplex::Vertex_handle const
+  >
+{
+	friend class boost::iterator_core_access;
+
+	typedef SkeletonBlockerComplex Complex;
+	typedef typename Complex::boost_adjacency_iterator boost_adjacency_iterator;
+	typedef typename Complex::Vertex_handle Vertex_handle;
+	typedef typename Complex::Edge_handle Edge_handle;
+
+private:
+
+	const Complex* complex;
+	Vertex_handle v;
+
+	boost_adjacency_iterator current_;
+	boost_adjacency_iterator end_;
+
+public:
+	//	boost_adjacency_iterator ai, ai_end;
+	//	for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end; ++ai){
+
+	Complex_neighbors_vertices_iterator():complex(NULL){
+	}
+
+	Complex_neighbors_vertices_iterator(const Complex* complex_,Vertex_handle v_):
+		complex(complex_),
+		v(v_){
+		tie(current_,end_) = adjacent_vertices(v.vertex, complex->skeleton);
+	}
+
+	/**
+	 * returns an iterator to the end
+	 */
+	Complex_neighbors_vertices_iterator(const Complex* complex_,Vertex_handle v_,int end):
+		complex(complex_),
+		v(v_){
+		tie(current_,end_) = adjacent_vertices(v.vertex, complex->skeleton);
+		set_end();
+	}
+
+
+	void increment () {
+		if(current_ != end_)
+			++current_;
+	}
+
+	Vertex_handle  dereference() const	{
+		return(Vertex_handle(*current_));
+	}
+
+	bool equal(const Complex_neighbors_vertices_iterator& other) const{
+		return (complex== other.complex) && (v == other.v) && (current_ == other.current_) && (end_ == other.end_);
+	}
+
+private:
+	//todo remove this ugly hack
+	void set_end(){
+		current_ = end_;
+	}
+};
+
+}
+
+} // namespace GUDHI
+
+#endif /* GUDHI_SKELETON_BLOCKERS_VERTICES_ITERATORS_H_ */
+
+
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
new file mode 100644
index 00000000..b4c1cd9e
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
@@ -0,0 +1,1422 @@
+/*
+ * Skeleton_blocker_complex.h
+ *
+ *  Created on: Jan, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKER_COMPLEX_H
+#define GUDHI_SKELETON_BLOCKER_COMPLEX_H
+
+
+#include <map>
+#include <list>
+#include <set>
+#include <vector>
+#include <iostream>
+#include <string>
+#include <fstream>
+#include <sstream>
+#include <memory>
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/connected_components.hpp>
+#include <boost/iterator/transform_iterator.hpp>
+#include <boost/range/adaptor/map.hpp>
+
+#include "gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h"
+#include "gudhi/Skeleton_blocker_link_complex.h"
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h"
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h"
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h"
+
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h"
+#include "gudhi/Skeleton_blocker/internal/Top_faces.h"
+#include "gudhi/Utils.h"
+
+
+namespace Gudhi {
+
+namespace skbl {
+
+
+/**
+ *@class Skeleton_blocker_complex
+ *@brief Abstract Simplicial Complex represented with a skeleton/blockers pair.
+ *
+ * A simplicial complex is completely defined by :
+ * - the graph of its 1-skeleton;
+ * - its set of blockers.
+ *
+ * The graph is a boost graph templated with SkeletonBlockerDS::Graph_vertex and SkeletonBlockerDS::Graph_edge.
+ *
+ * One can accesses to vertices via SkeletonBlockerDS::Vertex_handle, to edges via Skeleton_blocker_complex::Edge_handle and
+ * simplices via Skeleton_blocker_simplex<Vertex_handle>.
+ *
+ * The SkeletonBlockerDS::Root_vertex_handle serves in the case of a subcomplex (see class Skeleton_blocker_sub_complex)
+ * to access to the address of one vertex in the parent complex.
+ *
+ * @todo TODO Simplex_handle are not classic handle
+ *
+ */
+template<class SkeletonBlockerDS>
+class Skeleton_blocker_complex
+{
+	template<class ComplexType> friend class Complex_vertex_iterator;
+	template<class ComplexType> friend class Complex_neighbors_vertices_iterator;
+	template<class ComplexType> friend class Complex_edge_iterator;
+	template<class ComplexType> friend class Complex_edge_around_vertex_iterator;
+
+
+	template<class ComplexType> friend class Skeleton_blocker_link_complex;
+	template<class ComplexType> friend class Skeleton_blocker_link_superior;
+	template<class ComplexType> friend class Skeleton_blocker_sub_complex;
+
+public:
+
+
+	typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex;
+	typedef typename SkeletonBlockerDS::Graph_edge Graph_edge;
+
+	typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle;
+	typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle;
+	typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle;
+
+
+	typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle;
+	typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle;
+
+
+	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:
+	/**
+	 * Handle to an edge of the complex.
+	 */
+	typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle;
+
+
+
+	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. */
+
+
+	//todo remove!!!
+public:
+
+	/** Each vertex can access to the blockers passing through it. */
+	BlockerMap blocker_map_;
+
+
+
+
+public:
+
+
+
+
+	/////////////////////////////////////////////////////////////////////////////
+	/** @name Constructors / Destructors / Initialization
+	 */
+	//@{
+	Skeleton_blocker_complex(int num_vertices_ = 0,Visitor* visitor_=NULL):visitor(visitor_){
+		clear();
+		for (int i=0; i<num_vertices_; ++i){
+			add_vertex();
+		}
+	}
+
+private:
+
+	/**
+	 * 	this nested class is used for the next constructor that takes as an input a list of simplices.
+	 * 	It stores a vector where the ith cell contains a set of i-simplices
+	 *
+	 */
+	class Simplices_sets_from_list{
+	private:
+		typedef  std::set< Simplex_handle> Container_simplices;
+		typedef typename Container_simplices::iterator Simplices_iterator;
+		int dimension_;
+		std::vector<Container_simplices > simplices_;
+
+	public:
+		Simplices_sets_from_list(std::list<Simplex_handle>& simplices):
+			dimension_(-1){
+			assert(!simplices.empty());
+
+			for(auto simplex = simplices.begin() ; simplex != simplices.end(); ++simplex ){
+				dimension_ = std::max(dimension_,(int)simplex->dimension());
+			}
+			simplices_ = std::vector<Container_simplices >(dimension_+1);
+
+			// compute k-simplices
+			for(auto simplex = simplices.begin() ; simplex != simplices.end(); ++simplex ){
+				simplices_[simplex->dimension()].insert(*simplex);
+			}
+		}
+
+		Simplices_iterator begin(int k){
+			assert(0<= k && k<= dimension_);
+			return simplices_[k].begin();
+		}
+
+		Simplices_iterator end(int k){
+			assert(0<= k && k<= dimension_);
+			return simplices_[k].end();
+		}
+
+
+		Container_simplices& simplices(int k){
+			return simplices_[k];
+		}
+
+		int dimension(){
+			return dimension_;
+		}
+
+		bool contains(const Simplex_handle& simplex) const{
+			if(simplex.dimension()>dimension_)
+				return false;
+			else
+				return simplices_[simplex.dimension()].find(simplex)!= simplices_[simplex.dimension()].end();
+		}
+
+		void print(){
+			for(int i = 0; i < dimension_; ++i){
+				std::cout << i<<"-simplices"<<std::endl;
+				auto l = simplices_[i];
+				for(auto s : l){
+					std::cout << s<<std::endl;
+				}
+			}
+		}
+	};
+
+
+	void compute_next_expand(
+			Simplices_sets_from_list& simplices,
+			int dim,
+			std::list<Simplex_handle>& next_expand)
+	{
+		next_expand.clear();
+
+		//	    high_resolution_clock::time_point tbeginfor = high_resolution_clock::now();
+		//	    auto durationlooplink = std::chrono::duration_cast<std::chrono::microseconds>( tbeginfor - tbeginfor ).count();
+
+		for(auto sigma = simplices.begin(dim); sigma != simplices.end(dim); ++sigma){
+			//		    high_resolution_clock::time_point tbeg = high_resolution_clock::now();
+			Simplex_handle t(*sigma);
+			Skeleton_blocker_link_superior<Skeleton_blocker_complex> link(*this,t);
+			// xxx all time here, most likely because accessing superior edges needs passing through lower one
+			// currently
+
+			//			durationlooplink += std::chrono::duration_cast<std::chrono::microseconds>( high_resolution_clock::now() - tbeg ).count();
+
+			for(auto v : link.vertex_range()){
+				Vertex_handle v_in_complex(*this->get_address(  link.get_id(v)) );
+				t.add_vertex(v_in_complex);
+				next_expand.push_back(t);
+				t.remove_vertex(v_in_complex);
+			}
+		}
+		//	    high_resolution_clock::time_point t2 = high_resolution_clock::now();
+		//	    auto durationlooptotal = std::chrono::duration_cast<std::chrono::microseconds>( t2 - tbeginfor ).count();
+		//	    DBGVALUE(durationlooptotal);
+		//	    DBGVALUE(durationlooplink);
+	}
+
+
+
+public:
+	/**
+	 * @brief Constructor with a list of simplices
+	 * @details The list of simplices must be the list
+	 * of simplices of a simplicial complex.
+	 *todo rewrite as iterators
+	 */
+	Skeleton_blocker_complex(std::list<Simplex_handle>& simplices,Visitor* visitor_=NULL):
+		num_vertices_(0),num_blockers_(0),
+		visitor(visitor_){
+		Simplices_sets_from_list set_simplices(simplices);
+
+		int dim = set_simplices.dimension();
+
+
+		// add 1-skeleton to the complex
+		for(auto v_it = set_simplices.begin(0); v_it != set_simplices.end(0); ++v_it)
+			add_vertex();
+
+		for(auto e_it = set_simplices.begin(1); e_it != set_simplices.end(1); ++e_it){
+			Vertex_handle a = e_it->first_vertex();
+			Vertex_handle b = e_it->last_vertex();
+			assert(contains_vertex(a) && contains_vertex(b));
+			add_edge(a,b);
+		}
+
+		// then add blockers
+		for(int current_dim = 1 ; current_dim <=dim ; ++current_dim){
+			std::list<Simplex_handle> expansion_simplices;
+			compute_next_expand(set_simplices,current_dim,expansion_simplices);
+
+			for(const auto &simplex : expansion_simplices) {
+				if(!set_simplices.contains(simplex)){
+					add_blocker(simplex);
+				}
+			}
+		}
+	}
+
+
+
+	// 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;
+	}
+
+
+	/**
+	 * The destructor delete all blockers allocated.
+	 */
+	virtual ~Skeleton_blocker_complex(){
+		clear();
+	}
+
+	/**
+	 * 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();
+	}
+
+	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;
+	}
+
+	Graph_vertex& operator[](Vertex_handle address){
+		assert(0<=address.vertex && address.vertex< boost::num_vertices(skeleton));
+		return skeleton[address.vertex];
+	}
+
+	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 In fact, it just deactivates the vertex.
+	 */
+	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);
+	}
+
+	/**
+	 * @return true iff the simplicial complex contains the vertex u
+	 */
+	bool contains_vertex(Vertex_handle u) const{
+		if (u.vertex<0 || u.vertex>=boost::num_vertices(skeleton)) return false;
+		return (*this)[u].is_active();
+	}
+
+	/**
+	 * @return true iff the simplicial complex contains the vertex u
+	 */
+	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;
+	}
+
+	/**
+	 * Given an Id return the address of the vertex having this Id in the complex.
+	 * 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();
+	}
+
+
+	/**
+	 * 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
+	 */
+	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;
+	}
+
+
+	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;
+	}
+
+	Graph_edge& operator[](Edge_handle edge_handle){
+		return skeleton[edge_handle];
+	}
+
+	const Graph_edge& operator[](Edge_handle edge_handle) const{
+		return skeleton[edge_handle];
+	}
+
+	Vertex_handle first_vertex(Edge_handle edge_handle) const{
+		return source(edge_handle,skeleton);
+	}
+
+	Vertex_handle second_vertex(Edge_handle edge_handle) const{
+		return target(edge_handle,skeleton);
+	}
+
+	/**
+	 * @brief returns the simplex made with the two vertices of the edge
+	 */
+	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
+	 */
+	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 of simplex sigma 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 edge ab from the simplicial complex.
+	 */
+	virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b){
+		bool found;
+		Edge_handle edge;
+		tie(edge,found) = boost::edge(a.vertex,b.vertex,skeleton);
+		if (found)
+		{
+			if (visitor) visitor->on_remove_edge(a,b);
+			//		if (heapCollapse.Contains(edge)) heapCollapse.Delete(edge);
+			boost::remove_edge(a.vertex,b.vertex,skeleton);
+			degree_[a.vertex]--;
+			degree_[b.vertex]--;
+		}
+		return edge;
+	}
+
+
+	/**
+	 * @brief Removes edge 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
+	 */
+	//@{
+	/**
+	 * Adds the 2-blocker abc
+	 */
+	void add_blocker(Vertex_handle a, Vertex_handle b, Vertex_handle c){
+		add_blocker(Simplex_handle(a,b,c));
+	}
+
+	/**
+	 * Adds the 3-blocker abcd
+	 */
+	void add_blocker(Vertex_handle a, Vertex_handle b, Vertex_handle c, Vertex_handle d){
+		add_blocker(Simplex_handle(a,b,c,d));
+	}
+
+	/**
+	 * Adds the simplex blocker_pt to the set of blockers and
+	 * returns a Blocker_handle toward it if was not present before.
+	 */
+	Blocker_handle add_blocker(const Simplex_handle& blocker){
+		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:
+	/**
+	 * Adds the simplex s 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:
+	/**
+	 * Removes the simplex sigma from the set of blockers.
+	 * 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_--;
+	}
+
+
+
+
+	/**
+	 * 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 (auto blocker : const_blocker_range(a)){
+			if ( s == *blocker )
+				return true;
+		}
+		return false;
+	}
+
+	/**
+	 * @return true iff s is a blocker of the simplicial complex
+	 */
+	bool contains_blocker(const Simplex_handle & s) const{
+		if (s.dimension()<2)
+			return false;
+
+		Vertex_handle a = s.first_vertex();
+
+		for (auto blocker : const_blocker_range(a)){
+			if ( s == *blocker )
+				return true;
+		}
+		return false;
+	}
+
+
+private:
+	/**
+	 * @return true iff a blocker of the simplicial complex
+	 * is a face of sigma.
+	 */
+	bool blocks(const Simplex_handle & sigma) const{
+
+		for(auto blocker : const_blocker_range()){
+			if ( sigma.contains(*blocker) )
+				return true;
+		}
+		return false;
+	}
+
+	//@}
+
+
+
+
+	/////////////////////////////////////////////////////////////////////////////
+	/** @name Neighbourhood access
+	 */
+	//@{
+
+public:
+	/**
+	 * @brief Adds to simplex n the neighbours of v:
+	 * \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));
+		}
+	}
+
+	/**
+	 * @brief 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();
+		// ----------------------------
+		// Compute vertices in the link
+		// we compute the intersection of N(alpha_i) and store it in n
+		// ----------------------------
+		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);
+		}
+	}
+
+	/**
+	 * @brief Eliminates from simplex n all vertices which are
+	 * not neighbours of v: \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);
+	}
+
+	/**
+	 * @brief Eliminates from simplex n all vertices which are
+	 * neighbours of v: \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);
+	}
+	//@}
+
+
+	/////////////////////////////////////////////////////////////////////////////
+	/** @name Operations on the simplicial complex
+	 */
+	//@{
+public:
+
+	/**
+	 * @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.
+	 * todo in O(n), cache the value
+	 */
+	int num_edges() const{
+		return boost::num_edges(skeleton);
+	}
+
+	/*
+	 * @brief returns the number of blockers in the complex.
+	 */
+	int num_blockers() const{
+		return num_blockers_;
+	}
+
+	/*
+	 * @brief returns true iff the graph of the 1-skeleton of the complex is complete.
+	 */
+	bool complete() const{
+		return (num_vertices()*(num_vertices()-1))/2 == num_edges();
+	}
+
+	/**
+	 * @brief returns the number of connected components in the graph of the 1-skeleton.
+	 */
+	int num_connected_components() const{
+		int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices();
+		std::vector<int> component(skeleton.vertex_set().size());
+		return boost::connected_components(this->skeleton,&component[0]) - num_vert_collapsed;
+	}
+
+
+	//todo remove
+	// do
+	void keep_only_largest_cc(){
+		std::vector<unsigned> component(skeleton.vertex_set().size());
+		boost::connected_components(this->skeleton,&component[0]);
+		auto maxCC = min_element(component.begin(),component.end());
+		for (unsigned i = 0; i != component.size(); ++i){
+			if(component[i]!=*maxCC){
+				if(this->contains_vertex(Vertex_handle(i)))
+					this->remove_vertex(Vertex_handle(i));
+			}
+		}
+	}
+
+
+	/**
+	 * @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;
+	}
+
+	//@}
+
+
+	/////////////////////////////////////////////////////////////////////////////
+	/** @name Vertex iterators
+	 */
+	//@{
+
+	typedef Complex_vertex_iterator<Skeleton_blocker_complex> CVI; //todo rename
+
+	//	/**
+	//	 * @brief Range over the vertices of the simplicial complex.
+	//	 * Methods .begin() and .end() return a Complex_vertex_iterator.
+	//	 */
+	typedef boost::iterator_range     < Complex_vertex_iterator<Skeleton_blocker_complex> > Complex_vertex_range;
+
+	/**
+	 * @brief Returns a Complex_vertex_range over all vertices of the complex
+	 */
+	Complex_vertex_range vertex_range() const
+	{
+		auto begin = Complex_vertex_iterator<Skeleton_blocker_complex>(this);
+		auto end = Complex_vertex_iterator<Skeleton_blocker_complex>(this,0);
+		return Complex_vertex_range(begin,end);
+	}
+
+	typedef boost::iterator_range     < Complex_neighbors_vertices_iterator<Skeleton_blocker_complex> > Complex_neighbors_vertices_range;
+
+	/**
+	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v
+	 */
+	Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const
+	{
+		auto begin = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(this,v);
+		auto end = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(this,v,0);
+		return Complex_neighbors_vertices_range(begin,end);
+	}
+
+	//@}
+
+
+	/** @name Edge iterators
+	 */
+	//@{
+
+
+	typedef boost::iterator_range  <Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>>
+			Complex_edge_range;
+
+	/**
+	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex
+	 */
+	Complex_edge_range edge_range() const
+	{
+		auto begin = Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this);
+		auto end = Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this,0);
+		return Complex_edge_range(begin,end);
+	}
+
+
+
+	typedef boost::iterator_range  <Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>>
+			Complex_edge_around_vertex_range;
+	/**
+	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes
+	 * through 'v'
+	 */
+	Complex_edge_around_vertex_range edge_range(Vertex_handle v) const
+	{
+		auto begin = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this,v);
+		auto end = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this,v,0);
+		return Complex_edge_around_vertex_range(begin,end);
+	}
+
+
+
+	//@}
+
+
+	/** @name Triangles iterators
+	 */
+	//@{
+private:
+	typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link;
+	typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link;
+public:
+	typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex,Superior_link> Superior_triangle_around_vertex_iterator;
+
+	typedef boost::iterator_range     < Triangle_around_vertex_iterator<Skeleton_blocker_complex,Link> > Complex_triangle_around_vertex_range;
+
+	/**
+	 * @brief Range over triangles around a vertex of the simplicial complex.
+	 * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
+	 *
+	 */
+	Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const
+	{
+		auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex,Link>(this,v);
+		auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex,Link>(this,v,0);
+		return Complex_triangle_around_vertex_range(begin,end);
+	}
+
+
+	typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range;
+
+	/**
+	 * @brief Range over triangles of the simplicial complex.
+	 * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
+	 *
+	 */
+	Complex_triangle_range triangle_range() const
+	{
+		auto end = Triangle_iterator<Skeleton_blocker_complex>(this,0);
+		if(empty()){
+			return Complex_triangle_range(end,end);
+		}
+		else{
+			auto begin = Triangle_iterator<Skeleton_blocker_complex>(this);
+			return Complex_triangle_range(begin,end);
+		}
+
+	}
+
+
+	//@}
+
+
+
+	/** @name Simplices iterators
+	 */
+	//@{
+	typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex,Link> Complex_simplex_around_vertex_iterator;
+
+	/**
+	 * @brief Range over the simplices of the simplicial complex around a vertex.
+	 * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator.
+	 */
+	typedef boost::iterator_range     < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range;
+
+	/**
+	 * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex
+	 */
+	Complex_simplex_around_vertex_range simplex_range(Vertex_handle v) const
+	{
+		assert(contains_vertex(v));
+		return Complex_simplex_around_vertex_range(
+				Complex_simplex_around_vertex_iterator(this,v),
+				Complex_simplex_around_vertex_iterator(this,v,true)
+		);
+	}
+
+	//	typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator;
+	typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator;
+
+	typedef boost::iterator_range     < Complex_simplex_iterator > Complex_simplex_range;
+
+	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:
+
+
+
+	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 Complex_blocker_around_vertex_range over all blockers of the complex adjacent to the vertex v.
+	 */
+	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:
+
+
+	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);
+	}
+
+
+	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() << ")"<<" id = "<< (*this)[edge].index()<< std::endl;
+		stream<< std::endl;
+		return stream.str();
+	}
+
+
+	std::string blockers_to_string() const{
+		std::ostringstream stream;
+		for (auto bl:blocker_map_){
+			stream << bl.first << " => " << bl.second << ":"<<*bl.second <<"\n";
+		}
+		return stream.str();
+	}
+
+	//@}
+
+};
+
+
+
+/**
+ * build a simplicial complex from a collection
+ * of top faces.
+ */
+template<typename Complex,typename SimplexHandleIterator>
+unsigned make_complex_from_top_faces(Complex& complex,SimplexHandleIterator begin,SimplexHandleIterator end){
+	typedef typename Complex::Simplex_handle Simplex_handle;
+
+	int dimension = 0;
+	for(auto top_face = begin; top_face != end; ++top_face)
+		dimension = std::max(dimension,top_face->dimension());
+
+	std::vector< std::set<Simplex_handle> > simplices_per_dimension(dimension+1);
+
+
+	for(auto top_face = begin; top_face != end; ++top_face){
+		register_faces(simplices_per_dimension,*top_face);
+	}
+
+	// compute list of simplices
+	std::list<Simplex_handle> simplices;
+	for(int dim = 0 ; dim <= dimension ; ++dim){
+		std::copy(
+				simplices_per_dimension[dim].begin(),
+				simplices_per_dimension[dim].end(),
+				std::back_inserter(simplices)
+		);
+	}
+	complex = Complex(simplices);
+	return simplices.size();
+}
+
+
+
+}
+
+}  // namespace GUDHI
+
+
+
+#endif
+
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
new file mode 100644
index 00000000..924b653c
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
@@ -0,0 +1,118 @@
+/*
+ * Skeleton_blocker_geometric_complex.h
+ *
+ *  Created on: Feb 11, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_
+#define SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_
+
+
+#include "gudhi/Utils.h"
+#include "gudhi/Skeleton_blocker_simplifiable_complex.h"
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h"
+
+
+
+namespace Gudhi{
+
+namespace skbl {
+
+/**
+ * @brief Class that represents a geometric complex that can be simplified.
+ * The class allows access to points of vertices.
+ *
+ */
+template<typename SkeletonBlockerGeometricDS>
+class Skeleton_blocker_geometric_complex : public Skeleton_blocker_simplifiable_complex<SkeletonBlockerGeometricDS>
+{
+public:
+	typedef typename SkeletonBlockerGeometricDS::GT GT;
+
+	typedef Skeleton_blocker_simplifiable_complex<SkeletonBlockerGeometricDS> SimplifiableSkeletonblocker;
+
+	typedef typename SimplifiableSkeletonblocker::Vertex_handle Vertex_handle;
+	typedef typename SimplifiableSkeletonblocker::Root_vertex_handle Root_vertex_handle;
+	typedef typename SimplifiableSkeletonblocker::Edge_handle Edge_handle;
+	typedef typename SimplifiableSkeletonblocker::Simplex_handle Simplex_handle;
+
+
+	typedef typename SimplifiableSkeletonblocker::Graph_vertex Graph_vertex;
+
+	typedef typename SkeletonBlockerGeometricDS::Point Point;
+
+
+	/**
+	 * @brief Add a vertex to the complex with its associated point.
+	 */
+	void add_vertex(const Point& point){
+		Vertex_handle ad = SimplifiableSkeletonblocker::add_vertex();
+		(*this)[ad].point() = point;
+	}
+
+
+	/**
+	 * @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(const Simplex_handle& simplex) const{
+		Geometric_link link(*this,simplex);
+		//we now add the point info
+		add_points_to_link(link);
+		return link;
+	}
+
+	/**
+	 * Constructs the link of 'simplex' with points coordinates.
+	 */
+	Geometric_link link(Edge_handle edge) const{
+		Geometric_link link(*this,edge);
+		//we now add the point info
+		add_points_to_link(link);
+		return link;
+	}
+private:
+	void add_points_to_link(Geometric_link& link) const{
+		for(Vertex_handle v : link.vertex_range()){
+			Root_vertex_handle v_root(link.get_id(v));
+			link.point(v) = (*this).point(v_root);
+		}
+	}
+};
+
+}
+
+}  // namespace GUDHI
+
+#endif /* SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_ */
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_link_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_link_complex.h
new file mode 100644
index 00000000..29c0691b
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_link_complex.h
@@ -0,0 +1,275 @@
+/*
+ * Skeleton_blockers_link_complex.h
+ *
+ *  Created on: Feb 7, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKERS_LINK_COMPLEX_H_
+#define GUDHI_SKELETON_BLOCKERS_LINK_COMPLEX_H_
+
+#include "gudhi/Utils.h"
+#include "gudhi/Skeleton_blocker_complex.h"
+
+namespace Gudhi{
+
+namespace skbl {
+
+template<class ComplexType> class Skeleton_blocker_sub_complex;
+
+
+/**
+ *  \brief Class representing the link of a simplicial complex encoded by a skeleton/blockers pair.
+ *  It inherits from Skeleton_blocker_sub_complex because such complex is a sub complex of a
+ *  root complex.
+ */
+template<typename ComplexType>
+class Skeleton_blocker_link_complex : public Skeleton_blocker_sub_complex<ComplexType>
+{
+	template <typename T> friend class Skeleton_blocker_link_superior;
+	typedef typename ComplexType::Edge_handle Edge_handle;
+
+	typedef typename ComplexType::boost_vertex_handle boost_vertex_handle;
+
+
+private:
+
+	bool only_superior_vertices_;
+
+public:
+	typedef typename ComplexType::Vertex_handle Vertex_handle;
+	typedef typename ComplexType::Root_vertex_handle Root_vertex_handle;
+
+	typedef typename ComplexType::Simplex_handle Simplex_handle;
+	typedef typename ComplexType::Root_simplex_handle Root_simplex_handle;
+
+	typedef typename ComplexType::Blocker_handle Blocker_handle;
+
+
+	typedef typename ComplexType::Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator;
+	typedef typename ComplexType::Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_handle_iterator;
+
+
+	Skeleton_blocker_link_complex(bool only_superior_vertices=false):only_superior_vertices_(only_superior_vertices){
+	}
+
+	/**
+	 * If the parameter only_superior_vertices is true,
+	 * only vertices greater than the one of alpha are added.
+	 */
+	Skeleton_blocker_link_complex(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool only_superior_vertices = false)
+	:only_superior_vertices_(only_superior_vertices) {
+		if(!alpha_parent_adress.empty())
+			build_link(parent_complex,alpha_parent_adress);
+
+	}
+
+	/**
+	 * If the parameter only_superior_vertices is true,
+	 * only vertices greater than the one of the vertex are added.
+	 */
+	Skeleton_blocker_link_complex(const ComplexType & parent_complex, Vertex_handle a_parent_adress, bool only_superior_vertices = false)
+	:only_superior_vertices_(only_superior_vertices){
+		Simplex_handle alpha_simplex(a_parent_adress);
+		build_link(parent_complex,alpha_simplex);
+	}
+
+
+	/**
+	 * If the parameter only_superior_vertices is true,
+	 * only vertices greater than the one of the edge are added.
+	 */
+	Skeleton_blocker_link_complex(const ComplexType & parent_complex, Edge_handle edge, bool only_superior_vertices = false)
+	:only_superior_vertices_(only_superior_vertices){
+		Simplex_handle alpha_simplex(parent_complex.first_vertex(edge),parent_complex.second_vertex(edge));
+		build_link(parent_complex,alpha_simplex);
+	}
+
+protected:
+
+
+	/**
+	 * @brief compute vertices of the link.
+	 * If the boolean only_superior_vertices is true, then only the vertices
+	 * are greater than  vertices of alpha_parent_adress are added.
+	 */
+	void compute_link_vertices(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool only_superior_vertices,bool is_alpha_blocker = false){
+		if(alpha_parent_adress.dimension()==0)
+			// for a vertex we know exactly the number of vertices of the link (and the size of the corresponding vector)
+			// thus we call a specific function that will reserve a vector with appropriate size
+			this->compute_link_vertices(parent_complex,alpha_parent_adress.first_vertex(),only_superior_vertices_);
+		else{
+			// we compute the intersection of neighbors of alpha and store it in link_vertices
+			Simplex_handle link_vertices_parent;
+			parent_complex.add_neighbours(alpha_parent_adress,link_vertices_parent,only_superior_vertices);
+			// For all vertex 'v' in this intersection, we go through all its adjacent blockers.
+			// If one blocker minus 'v' is included in alpha then the vertex is not in the link complex.
+			for (auto v_parent : link_vertices_parent){
+				bool new_vertex = true;
+				for (auto beta : parent_complex.const_blocker_range(v_parent))
+				{
+					if(!is_alpha_blocker || *beta!=alpha_parent_adress){
+						new_vertex = !(alpha_parent_adress.contains_difference(*beta,v_parent));
+						if(!new_vertex) break;
+					}
+				}
+				if (new_vertex)
+					this->add_vertex(parent_complex.get_id(v_parent));
+			}
+		}
+	}
+
+
+	/**
+	 * @brief compute vertices of the link.
+	 * If the boolean only_superior_vertices is true, then only the vertices
+	 * are greater than  vertices of alpha_parent_adress are added.
+	 */
+	void compute_link_vertices(const ComplexType & parent_complex, Vertex_handle alpha_parent_adress,bool only_superior_vertices){
+		// for a vertex we know exactly the number of vertices of the link (and the size of the corresponding vector
+		this->skeleton.m_vertices.reserve(parent_complex.degree(alpha_parent_adress));
+
+		// For all vertex 'v' in this intersection, we go through all its adjacent blockers.
+		// If one blocker minus 'v' is included in alpha then the vertex is not in the link complex.
+		for (auto v_parent : parent_complex.vertex_range(alpha_parent_adress)){
+			if(!only_superior_vertices || v_parent.vertex> alpha_parent_adress.vertex)
+				this->add_vertex(parent_complex.get_id(v_parent));
+		}
+
+	}
+
+
+	void compute_link_edges(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool is_alpha_blocker = false){
+		Simplex_handle_iterator y_link, x_parent , y_parent;
+		// ----------------------------
+		// Compute edges in the link
+		// -------------------------
+
+		if(this->num_vertices()<=1) return;
+
+		for (auto x_link = this->vertex_range().begin() ;
+				x_link!= this->vertex_range().end();
+				++x_link
+		)
+		{
+			for (auto y_link = x_link; ++y_link != this->vertex_range().end(); ){
+				Vertex_handle x_parent = *parent_complex.get_address(this->get_id(*x_link));
+				Vertex_handle y_parent = *parent_complex.get_address(this->get_id(*y_link));
+				if (parent_complex.contains_edge(x_parent,y_parent) ){
+					// we check that there is no blocker subset of alpha passing trough x and y
+					bool new_edge=true;
+					for (auto blocker_parent : parent_complex.const_blocker_range(x_parent))
+					{
+						if(!is_alpha_blocker || *blocker_parent!=alpha_parent_adress){
+							if 	(blocker_parent->contains(y_parent))
+							{
+								new_edge = !(alpha_parent_adress.contains_difference(*blocker_parent,x_parent,y_parent));
+								if (!new_edge) break;
+							}
+						}
+					}
+					if (new_edge)
+						this->add_edge(*x_link,*y_link);
+				}
+			}
+		}
+	}
+
+
+	/**
+	 * @brief : Given an address in the current complex, it returns the
+	 * corresponding address in 'other_complex'.
+	 * It assumes that other_complex have a vertex 'this.get_id(address)'
+	 */
+	boost::optional<Vertex_handle> give_equivalent_vertex(const ComplexType & other_complex,
+			Vertex_handle address) const{
+		Root_vertex_handle id((*this)[address].get_id());
+		return other_complex.get_address(id);
+	}
+
+	/*
+	 * compute the blockers of the link if is_alpha_blocker is false.
+	 * Otherwise, alpha is a blocker, and the link is computed in the complex where
+	 * the blocker is anticollapsed.
+	 */
+	void compute_link_blockers(const ComplexType & parent_complex,const Simplex_handle& alpha_parent,bool is_alpha_blocker = false){
+
+		for (auto x_link : this->vertex_range()){
+
+			Vertex_handle x_parent = * this->give_equivalent_vertex(parent_complex,x_link);
+
+			for (auto blocker_parent : parent_complex.const_blocker_range(x_parent)){
+				if(!is_alpha_blocker || *blocker_parent!=alpha_parent){
+					Simplex_handle sigma_parent(*blocker_parent);
+
+					sigma_parent.difference(alpha_parent);
+
+					if (sigma_parent.dimension()>=2 && sigma_parent.first_vertex() == x_parent){
+						Root_simplex_handle sigma_id(parent_complex.get_id(sigma_parent));
+						auto  sigma_link = this->get_simplex_address(sigma_id);
+						// ie if the vertices of sigma are vertices of the link
+						if(sigma_link){
+							bool is_new_blocker = true;
+							for(auto a : alpha_parent){
+								for(auto eta_parent : parent_complex.const_blocker_range(a)){
+									if(!is_alpha_blocker || *eta_parent!=alpha_parent){
+										Simplex_handle eta_minus_alpha(*eta_parent);
+										eta_minus_alpha.difference(alpha_parent);
+										if(eta_minus_alpha != sigma_parent &&
+												sigma_parent.contains_difference(*eta_parent,alpha_parent)){
+											is_new_blocker = false;
+											break;
+										}
+									}
+								}
+								if (!is_new_blocker) break;
+							}
+							if (is_new_blocker)
+								this->add_blocker(new Simplex_handle(*sigma_link));
+
+						}
+					}
+				}
+			}
+		}
+	}
+
+
+public:
+	/**
+	 * @brief compute vertices, edges and blockers of the link.
+	 * @details If the boolean only_superior_vertices is true, then the link is computed only
+	 * with vertices that are greater than  vertices of alpha_parent_adress.
+	 */
+	void build_link(const ComplexType & parent_complex,const Simplex_handle& alpha_parent_adress,bool is_alpha_blocker = false)
+	{
+		assert(is_alpha_blocker || parent_complex.contains(alpha_parent_adress));
+		compute_link_vertices(parent_complex,alpha_parent_adress,only_superior_vertices_);
+		compute_link_edges(parent_complex,alpha_parent_adress,is_alpha_blocker);
+		compute_link_blockers(parent_complex,alpha_parent_adress,is_alpha_blocker);
+	}
+
+
+	/**
+	 * @brief build the link of a blocker which is the link
+	 * of the blocker's simplex if this simplex had been
+	 * removed from the blockers of the complex.
+	 */
+	friend void build_link_of_blocker(
+			const ComplexType & parent_complex,
+			Simplex_handle& blocker,
+			Skeleton_blocker_link_complex & result){
+		assert(blocker.dimension()>=2);
+		assert(parent_complex.contains_blocker(blocker));
+		result.clear();
+		result.build_link(parent_complex,blocker,true);
+	}
+
+
+};
+
+}
+
+}  // namespace GUDHI
+
+#endif /* SKELETON_BLOCKERS_LINK_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
new file mode 100644
index 00000000..4f9060ef
--- /dev/null
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
@@ -0,0 +1,525 @@
+/*
+ * Skeleton_blocker_simplifiable_complex.h
+ *
+ *  Created on: Feb 4, 2014
+ *      Author: dsalinas
+ */
+
+#ifndef GUDHI_SKELETON_BLOCKERS_SIMPLIFIABLE_COMPLEX_H_
+#define GUDHI_SKELETON_BLOCKERS_SIMPLIFIABLE_COMPLEX_H_
+
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h"
+
+namespace Gudhi{
+
+namespace skbl {
+
+/**
+ *  \brief Class that allows simplification operation on a simplicial complex represented
+ *  by a skeleton/blockers pair.
+ */
+template<typename SkeletonBlockerDS>
+class Skeleton_blocker_simplifiable_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
+	 */
+	//@{
+	Skeleton_blocker_simplifiable_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.
+	 * todo take iterator instead
+	 */
+	Skeleton_blocker_simplifiable_complex(std::list<Simplex_handle>& simplices,Visitor* visitor_=NULL):
+		Skeleton_blocker_complex<SkeletonBlockerDS>(simplices,visitor_)
+		{}
+
+
+
+	virtual ~Skeleton_blocker_simplifiable_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.
+	 *
+	 */
+	virtual bool is_popable_blocker(Blocker_handle sigma) const{
+		assert(this->contains_blocker(*sigma));
+		Skeleton_blocker_link_complex<Skeleton_blocker_simplifiable_complex> link_blocker_sigma;
+		build_link_of_blocker(*this,*sigma,link_blocker_sigma);
+
+		bool res = link_blocker_sigma.is_contractible()==CONTRACTIBLE;
+		return res;
+	}
+
+
+
+
+private:
+	/**
+	 * @returns the list of blockers of the simplex
+	 *
+	 * @todo a enlever et faire un iterateur sur tous les blockers a la place
+	 */
+	std::list<Blocker_handle> get_blockers(){
+		std::list<Blocker_handle> res;
+		for (auto blocker : this->blocker_range()){
+			res.push_back(blocker);
+		}
+		return res;
+	}
+
+
+
+public:
+
+	/**
+	 * Removes all the popable blockers of the complex and delete them.
+	 * @returns the number of popable blockers deleted
+	 */
+	void 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;
+					}
+				}
+			}
+		}
+	}
+
+
+	/**
+	 * Removes all the popable blockers of the complex passing through v and delete them.
+	 */
+	void 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;
+				}
+			}
+		}
+	}
+
+
+	/**
+	 * Remove the star of the vertex 'v'
+	 */
+	void remove_star(Vertex_handle v){
+		// we remove the blockers that are not consistent anymore
+
+		update_blockers_after_remove_star_of_vertex_or_edge(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);
+	}
+
+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){
+		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);
+		}
+	}
+
+
+
+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){
+		update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(a,b));
+		// we remove the edge
+		this->remove_edge(a,b);
+	}
+
+
+	/**
+	 * Remove the star of the edge 'e'.
+	 */
+	void remove_star(Edge_handle 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
+	 */
+	void 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
+				update_blockers_after_remove_star_of_simplex(sigma);
+	}
+
+private:
+	void update_blockers_after_remove_star_of_simplex(const Simplex_handle& sigma){
+		std::list <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);
+		this->add_blocker(sigma);
+	}
+
+public:
+
+	enum simplifiable_status{ NOT_HOMOTOPY_EQ,MAYBE_HOMOTOPY_EQ,HOMOTOPY_EQ};
+	simplifiable_status is_remove_star_homotopy_preserving(const Simplex_handle& simplex){
+		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;
+		//		return this->is_cone();
+	}
+
+
+	/** @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.
+	 * Suppose that the link condition Link(ab) = Link(a) inter Link(b)
+	 * is satisfied.
+	 */
+	void 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);
+		}
+	}
+
+
+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){
+		this->add_edge(a,x);
+		if (this->visitor) this->visitor->on_swaped_edge(a,b,x);
+		this->remove_edge(b,x);
+	}
+
+
+private:
+	/**
+	 * @brief removes all blockers passing through the edge 'ab'
+	 */
+	void 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();
+		}
+
+	}
+	/**
+	 * @brief removes all blockers passing through the edge 'ab'
+	 */
+	void 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();
+		}
+	}
+
+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){
+		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);
+
+		std::set<Simplex_handle> blockers_to_add;
+
+		get_blockers_to_be_added_after_contraction(a,b,blockers_to_add);
+
+		delete_blockers_around_vertices(a,b);
+
+		update_edges_after_contraction(a,b);
+
+		this->remove_vertex(b);
+
+		notify_changed_edges(a);
+
+		for(auto block : blockers_to_add)
+			this->add_blocker(block);
+
+		assert(this->contains_vertex(a));
+		assert(!this->contains_vertex(b));
+	}
+
+
+private:
+
+	void 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<SkeletonBlockerComplex>(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
+	 */
+	void 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();
+		}
+	}
+
+
+	void 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);
+		}
+	}
+
+	void 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));
+	}
+
+	//@}
+
+
+};
+
+}
+
+}  // namespace GUDHI
+
+#endif /* GUDHI_SKELETON_BLOCKERS_SIMPLIFIABLE_COMPLEX_H_ */
diff --git a/src/Skeleton_blocker/test/CMakeLists.txt b/src/Skeleton_blocker/test/CMakeLists.txt
new file mode 100644
index 00000000..d66dcf64
--- /dev/null
+++ b/src/Skeleton_blocker/test/CMakeLists.txt
@@ -0,0 +1,7 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIskbl)
+
+add_executable ( TestSkeletonBlockerComplex TestSkeletonBlockerComplex.cpp )
+add_executable ( TestSimplifiable           TestSimplifiable.cpp )
+
+
diff --git a/src/Skeleton_blocker/test/TestSimplifiable.cpp b/src/Skeleton_blocker/test/TestSimplifiable.cpp
new file mode 100644
index 00000000..9f802463
--- /dev/null
+++ b/src/Skeleton_blocker/test/TestSimplifiable.cpp
@@ -0,0 +1,281 @@
+/*
+ * TestSimplifiable.cxx
+ *
+ *  Created on: Feb 4, 2014
+ *      Author: dsalinas
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+#include <fstream>
+#include <sstream>
+#include "gudhi/Test.h"
+//#include "Skeleton_blocker/Simplex.h"
+#include "gudhi/Skeleton_blocker_complex.h"
+#include "gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h"
+#include "gudhi/Skeleton_blocker_simplifiable_complex.h"
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h"
+
+
+using namespace std;
+
+using namespace Gudhi;
+
+using namespace skbl;
+
+template<typename ComplexType> class Skeleton_blocker_sub_complex;
+typedef Skeleton_blocker_simplifiable_complex<Skeleton_blocker_simple_traits> Complex;
+typedef Complex::Vertex_handle Vertex_handle;
+typedef Complex::Root_vertex_handle Root_vertex_handle;
+typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle;
+// true iff v \in complex
+bool assert_vertex(Complex &complex,Vertex_handle v){
+	Simplex_handle simplex(v);
+	bool test = complex.contains(simplex);
+	assert(test);
+	return test;
+}
+
+// true iff the blocker (a,b,c) is in complex
+bool assert_blocker(Complex &complex,Root_vertex_handle a,Root_vertex_handle b,Root_vertex_handle c){
+
+	return complex.contains_blocker(Simplex_handle(*complex.get_address(a),*complex.get_address(b),*complex.get_address(c)));
+	//return complex.contains_blocker((a),(b),(c));
+}
+
+void build_complete(int n,Complex& complex){
+	complex.clear();
+	for(int i=0;i<n;i++)
+		complex.add_vertex();
+	for(int i=0;i<n;i++)
+		for(int j=0;j<i;j++)
+			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+}
+
+bool test_contraction1(){
+	enum { a, b, x, y, z, n };
+	Complex complex(n);
+	build_complete(n,complex);
+	complex.remove_edge(Vertex_handle(b),Vertex_handle(z));
+	complex.add_blocker(Vertex_handle(a),Vertex_handle(x),Vertex_handle(y));
+	complex.add_blocker(Vertex_handle(b),Vertex_handle(x),Vertex_handle(y));
+
+	// Print result
+	cerr << "complex before complex"<< complex.to_string()<<endl;
+
+	cerr <<endl<<endl;
+	complex.contract_edge(Vertex_handle(a),Vertex_handle(b));
+	// Print result
+	cerr << "ContractEdge(0,1)\n";
+	PRINT(complex.to_string());
+
+	// verification
+	for (int i=0;i<5;i++)
+		if (i!=1) assert_vertex(complex,Vertex_handle(i));
+	bool test1 = !complex.contains_edge(Vertex_handle(a),Vertex_handle(b));
+	bool test2 = assert_blocker(complex,Root_vertex_handle(a),Root_vertex_handle(x),Root_vertex_handle(y));
+	bool test3 = complex.num_edges()==6;
+	bool test4 = complex.num_blockers()==1;
+	Simplex_handle sigma;
+	sigma.add_vertex(Vertex_handle(a));
+	sigma.add_vertex(Vertex_handle(x));
+	sigma.add_vertex(Vertex_handle(y));
+	sigma.add_vertex(Vertex_handle(z));
+	bool test5 = !(complex.contains(sigma));
+	return test1&&test2&&test3&&test4&&test5;
+}
+
+
+bool test_contraction2(){
+	enum { a, b, x, y, z, n };
+	Complex complex(n);
+	build_complete(n,complex);
+	complex.remove_edge(Vertex_handle(b),Vertex_handle(x));
+	Simplex_handle blocker;
+	blocker.add_vertex(Vertex_handle(a));
+	blocker.add_vertex(Vertex_handle(y));
+	blocker.add_vertex(Vertex_handle(z));
+
+	complex.add_blocker(blocker);
+
+	// Print result
+	cerr << "complex complex"<< complex.to_string();
+	cerr <<endl<<endl;
+	complex.contract_edge(Vertex_handle(a),Vertex_handle(b));
+
+	cerr << "complex.ContractEdge(a,b)"<< complex.to_string();
+
+	cerr <<endl<<endl;
+
+	// there should be one blocker (a,c,d,e) in the complex
+	bool test ;
+	test = complex.contains_blocker(Simplex_handle(a,x,y,z));
+	test = test && complex.num_blockers()==1;
+	return test;
+}
+
+bool test_link_condition1(){
+
+	Complex complex(0);
+	// Build the complexes
+	build_complete(4,complex);
+	complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2));
+
+
+	// Print result
+	cerr << "complex complex"<< complex.to_string();
+	cerr <<endl<<endl;
+
+	bool weak_link_condition = complex.link_condition(Vertex_handle(1),Vertex_handle(2),true);
+
+	bool strong_link_condition = complex.link_condition(Vertex_handle(1),Vertex_handle(2),false);
+
+	return weak_link_condition && !strong_link_condition;
+}
+
+
+
+
+bool test_collapse1(){
+	// xxx implement remove_star(simplex) before
+
+	Complex complex(5);
+	build_complete(4,complex);
+	complex.add_vertex();
+	complex.add_edge(2,4);
+	complex.add_edge(3,4);
+	// Print result
+	cerr << "initial complex :\n"<< complex.to_string();
+	cerr <<endl<<endl;
+
+	Simplex_handle simplex_123(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3));
+	complex.remove_star(simplex_123);
+	cerr << "complex.remove_star(1,2,3):\n"<< complex.to_string();
+	cerr <<endl<<endl;
+
+	// verification
+	bool blocker123_here = complex.contains_blocker(simplex_123);
+	cerr <<"----> Ocomplex \n";
+	return blocker123_here;
+}
+
+bool test_collapse2(){
+	Complex complex(5);
+	build_complete(4,complex);
+	complex.add_vertex();
+	complex.add_edge(Vertex_handle(1),Vertex_handle(4));
+	complex.add_edge(Vertex_handle(2),Vertex_handle(4));
+	complex.add_edge(Vertex_handle(3),Vertex_handle(4));
+	complex.add_blocker(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3),Vertex_handle(4));
+	// Print result
+	cerr << "initial complex :\n"<< complex.to_string();
+	cerr <<endl<<endl;
+
+	Simplex_handle sigma(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3));
+	complex.remove_star(sigma);
+	cerr << "complex.remove_star(1,2,3):\n"<< complex.to_string();
+	cerr <<endl<<endl;
+
+	// verification
+	bool blocker_removed = !complex.contains_blocker(Simplex_handle(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3),Vertex_handle(4)));
+	bool blocker123_here = complex.contains_blocker(sigma);
+	return blocker_removed && blocker123_here;
+}
+
+bool test_collapse3(){
+	Complex complex(5);
+	build_complete(4,complex);
+	complex.add_vertex();
+	complex.add_edge(Vertex_handle(1),Vertex_handle(4));
+	complex.add_edge(Vertex_handle(2),Vertex_handle(4));
+	complex.add_edge(Vertex_handle(3),Vertex_handle(4));
+	complex.add_blocker(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3),Vertex_handle(4));
+	// Print result
+	cerr << "initial complex:\n"<< complex.to_string();
+	cerr <<endl<<endl;
+
+	complex.remove_star(Vertex_handle(2));
+	cerr << "complex after remove star of 2:\n"<< complex.to_string();
+
+	bool blocker134_here = complex.contains_blocker(Simplex_handle(Vertex_handle(1),Vertex_handle(3),Vertex_handle(4)));
+	bool blocker1234_here = complex.contains_blocker(Simplex_handle(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3),Vertex_handle(4)));
+	// verification
+	//	assert_blocker(complex,1,2,3);
+	//	assert(!complex.ContainsBlocker(new AddressSimplex(1,2,3,4)));
+	cerr <<"----> Ocomplex \n";
+	return blocker134_here && !blocker1234_here;
+}
+
+
+
+bool test_remove_popable_blockers(){
+	Complex complex;
+	build_complete(4,complex);
+	complex.add_vertex();
+	complex.add_edge(Vertex_handle(3),Vertex_handle(4));
+	complex.add_edge(Vertex_handle(2),Vertex_handle(4));
+	Simplex_handle sigma1=Simplex_handle(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3));
+	Simplex_handle sigma2=Simplex_handle(Vertex_handle(2),Vertex_handle(3),Vertex_handle(4));
+
+	complex.add_blocker(sigma1);
+	complex.add_blocker(sigma2);
+	cerr << "complex complex"<< complex.to_string();
+	cerr <<endl<<endl;
+	cerr << "complex.RemovePopableBlockers();"<<endl;
+	complex.remove_popable_blockers();
+	cerr << "complex complex"<< complex.to_string();
+	cerr <<endl<<endl;
+
+	bool test1 = (complex.num_blockers()==1);
+
+
+	// test 2
+	complex.clear();
+	build_complete(4,complex);
+	complex.add_vertex();
+	complex.add_vertex();
+	complex.add_edge(Vertex_handle(3),Vertex_handle(4));
+	complex.add_edge(Vertex_handle(2),Vertex_handle(4));
+	complex.add_edge(Vertex_handle(2),Vertex_handle(5));
+	complex.add_edge(Vertex_handle(3),Vertex_handle(5));
+	complex.add_edge(Vertex_handle(4),Vertex_handle(5));
+	sigma1=Simplex_handle(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3));
+	sigma2=Simplex_handle(Vertex_handle(2),Vertex_handle(3),Vertex_handle(4));
+
+	complex.add_blocker(sigma1);
+	complex.add_blocker(sigma2);
+	cerr << "complex complex"<< complex.to_string();
+	cerr <<endl<<endl;	cerr << "complex.RemovePopableBlockers();"<<endl;
+	complex.remove_popable_blockers();
+	cerr << "complex complex"<< complex.to_string();
+
+	cerr <<endl<<endl;
+	bool test2 = (complex.num_blockers()==0);
+	return test1&&test2;
+}
+
+
+
+int main (int argc, char *argv[])
+{
+	Tests tests_simplifiable_complex;
+	tests_simplifiable_complex.add("Test contraction 1",test_contraction1);
+	tests_simplifiable_complex.add("Test contraction 2",test_contraction2);
+	tests_simplifiable_complex.add("Test Link condition 1",test_link_condition1);
+	tests_simplifiable_complex.add("Test remove popable blockers",test_remove_popable_blockers);
+
+
+	tests_simplifiable_complex.add("Test collapse 1",test_collapse1);
+	tests_simplifiable_complex.add("Test collapse 2",test_collapse2);
+	tests_simplifiable_complex.add("Test collapse 3",test_collapse3);
+
+
+	tests_simplifiable_complex.run();
+
+	if(tests_simplifiable_complex.run())
+		return EXIT_SUCCESS;
+	else
+		return EXIT_FAILURE;
+}
diff --git a/src/Skeleton_blocker/test/TestSkeletonBlockerComplex.cpp b/src/Skeleton_blocker/test/TestSkeletonBlockerComplex.cpp
new file mode 100644
index 00000000..c89f84e2
--- /dev/null
+++ b/src/Skeleton_blocker/test/TestSkeletonBlockerComplex.cpp
@@ -0,0 +1,773 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+#include <fstream>
+#include <sstream>
+#include "gudhi/Utils.h"
+#include "gudhi/Test.h"
+//#include "Skeleton_blocker/Simplex.h"
+#include "gudhi/Skeleton_blocker_complex.h"
+#include "gudhi/Skeleton_blocker_link_complex.h"
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h"
+#include "gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h"
+//#include "Simple_vertex.h"
+//#include "Simple_edge.h"
+
+using namespace std;
+
+using namespace Gudhi;
+
+using namespace skbl;
+
+
+typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> Complex;
+typedef Complex::Vertex_handle Vertex_handle;
+typedef Complex::Root_vertex_handle Root_vertex_handle;
+typedef Complex::Simplex_handle Simplex_handle;
+typedef Complex::Root_simplex_handle Root_simplex_handle;
+typedef Simplex_handle::Simplex_vertex_const_iterator Simplex_vertex_const_iterator;
+typedef Complex::Edge_handle Edge_handle;
+
+// true iff v \in complex
+bool assert_vertex(Complex &complex,Vertex_handle v){
+	//assert(complex.contains(v));
+	return complex.contains(v);
+}
+
+bool assert_simplex(Complex &complex,Root_vertex_handle a,Root_vertex_handle b,Root_vertex_handle c){
+	return true;
+	//	AddressSimplex simplex((a),(b),(c));
+	//	return complex.contains(&simplex);
+}
+
+// true iff the blocker (a,b,c) is in complex
+bool assert_blocker(Complex &complex,Root_vertex_handle a,Root_vertex_handle b,Root_vertex_handle c){
+	return true;
+	//return complex.contains_blocker((a),(b),(c));
+}
+
+// true iff the blocker (a,b,c,d) is in complex
+bool assert_blocker(Complex &complex,Root_vertex_handle a,Root_vertex_handle b,Root_vertex_handle c,Root_vertex_handle d){
+	return true;
+	//Simplex blocker (a,b,c,d);
+	//return complex.contains_blocker(&blocker);
+}
+
+
+void build_complete(int n,Complex& complex){
+	complex.clear();
+	for(int i=0;i<n;i++)
+		complex.add_vertex();
+
+//	for(int i=n-1;i>=0;i--)
+//		for(int j=i-1;j>=0;j--)
+//			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+
+	for(int i=0;i<n;i++)
+		for(int j=0;j<i;j++)
+			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+}
+
+
+bool test_simplex(){
+//	PRINT("test simplex");
+	Simplex_handle simplex(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2),Vertex_handle(3));
+	for (auto i = simplex.begin() ; i != simplex.end() ; ++i){
+		PRINT(*i);
+		auto j = i;
+		for (++j ;
+				j != simplex.end() ;
+				++j){
+			PRINT(*j);
+		}
+	}
+	return simplex.dimension()==3;
+}
+
+
+bool test_iterator_vertices1(){
+	int n = 10;
+	Complex complex(10);
+	cerr << "complex.num_vertices():"<<complex.num_vertices()<<endl;
+	int num_vertex_seen = 0;
+	for(auto vi :complex.vertex_range()){
+		cerr << "vertex:"<<vi<<endl;
+		++num_vertex_seen;
+	}
+	return num_vertex_seen == n;
+}
+
+bool test_iterator_vertices2(){
+	int n = 10;
+	Complex complex(10);
+	build_complete(10,complex);
+	cerr << "complex.num_vertices():"<<complex.num_vertices()<<endl;
+	cerr << "complex.num_edges():"<<complex.num_edges()<<endl;
+	int num_vertex_seen = 0;
+	for(auto vi :complex.vertex_range(Vertex_handle(2))){
+		cerr << "vertex:"<<vi<<endl;
+		++num_vertex_seen;
+	}
+	std::cerr<<"num_vertex_seen:"<<num_vertex_seen<<std::endl;
+	return num_vertex_seen == (n-1);
+}
+
+
+
+bool test_iterator_edge(){
+	const int n = 10;
+	Complex complex(n);
+	for(int i=0;i<n;i++)
+		for(int j=0;j<i;j++)
+			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+	complex.remove_edge(Vertex_handle(2),Vertex_handle(3));
+	complex.remove_edge(Vertex_handle(3),Vertex_handle(5));
+	cerr << "complex.num_edges():"<<complex.num_edges()<<endl;
+	int num_edges_seen = 0;
+	for(auto edge : complex.edge_range()){
+		cerr << "edge :"<<complex[edge]<<endl;
+		++num_edges_seen;
+	}
+
+	return num_edges_seen == n*(n-1)/2-2;
+}
+
+bool test_iterator_edge2(){
+	const int n = 10;
+	Complex complex(n);
+	for(int i=0;i<n;i++)
+		for(int j=0;j<i;j++)
+			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+	complex.remove_edge(Vertex_handle(2),Vertex_handle(3));
+	complex.remove_edge(Vertex_handle(3),Vertex_handle(5));
+	cerr << "complex.num_edges():"<<complex.num_edges()<<endl;
+	int num_neigbors_seen = 0;
+	for(auto neighbor : complex.vertex_range(Vertex_handle(2))){
+		cerr << "neighbor"<<neighbor<<endl;
+		++num_neigbors_seen;
+	}
+	return num_neigbors_seen==8;
+}
+
+
+
+bool test_iterator_edge3(){
+	const int n = 10;
+	Complex complex(n);
+	for(int i=0;i<n;i++)
+		for(int j=0;j<i;j++)
+			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+	complex.remove_edge(Vertex_handle(2),Vertex_handle(3));
+	complex.remove_edge(Vertex_handle(3),Vertex_handle(5));
+	cerr << "complex.num_edges():"<<complex.num_edges()<<endl;
+	int num_neigbors_seen = 0;
+	for(auto edge : complex.edge_range(Vertex_handle(2))){
+		std::cerr << edge<< std::endl;
+		++num_neigbors_seen;
+	}
+	return num_neigbors_seen==8;
+}
+
+
+
+bool test_iterator_triangles(){
+	const int n = 7;
+	Complex complex(n);
+	//create a "ring" around '0'
+	for(int i=1;i<n;i++)
+		complex.add_edge(Vertex_handle(0),Vertex_handle(i));
+	for(int i=1;i<n-1;i++)
+		complex.add_edge(Vertex_handle(i),Vertex_handle(i+1));
+	complex.add_edge(Vertex_handle(1),Vertex_handle(6));
+
+	PRINT(complex.to_string());
+
+	int num_triangles_seen=0;
+	//for (auto t : complex.triangle_range(5)){
+	TEST("triangles around 5 (should be 2 of them):");
+	for (auto t : complex.triangle_range(5)){
+		PRINT(t);
+		++num_triangles_seen;
+	}
+	bool test = (num_triangles_seen==2);
+
+	num_triangles_seen=0;
+	TEST("triangles around 0 (should be 6 of them):");
+	for (auto t : complex.triangle_range(0)){
+		PRINT(t);
+		++num_triangles_seen;
+	}
+	test = test&&(num_triangles_seen==6);
+
+	// we now add another triangle
+	complex.add_vertex();
+	complex.add_edge(Vertex_handle(4),Vertex_handle(7));
+	complex.add_edge(Vertex_handle(3),Vertex_handle(7));
+	complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(6));
+	num_triangles_seen=0;
+
+	TEST("triangles (should be 6 of them):");
+	num_triangles_seen=0;
+	for (auto t : complex.triangle_range()){
+		PRINT(t);
+		++num_triangles_seen;
+	}
+	test = test&&(num_triangles_seen==6);
+	PRINT(num_triangles_seen);
+
+	return test;
+}
+
+
+//#include "combinatorics/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h"
+
+bool test_iterator_simplices(){
+	Complex complex(6);
+	complex.add_edge(Vertex_handle(0),Vertex_handle(1));
+	complex.add_edge(Vertex_handle(1),Vertex_handle(2));
+	complex.add_edge(Vertex_handle(2),Vertex_handle(0));
+	complex.add_edge(Vertex_handle(1),Vertex_handle(3));
+	complex.add_edge(Vertex_handle(2),Vertex_handle(3));
+	complex.add_edge(Vertex_handle(2),Vertex_handle(5));
+	complex.add_edge(Vertex_handle(3),Vertex_handle(5));
+	complex.add_edge(Vertex_handle(2),Vertex_handle(4));
+	complex.add_edge(Vertex_handle(4),Vertex_handle(5));
+	complex.add_edge(Vertex_handle(3),Vertex_handle(4));
+
+	complex.add_blocker(Vertex_handle(2),Vertex_handle(3),Vertex_handle(4),Vertex_handle(5));
+
+	bool correct_number_simplices = true;
+
+	std::map<Vertex_handle,unsigned> expected_num_simplices;
+
+	expected_num_simplices[Vertex_handle(0)] = 4;
+	expected_num_simplices[Vertex_handle(1)] = 6;
+	expected_num_simplices[Vertex_handle(2)] = 11;
+	expected_num_simplices[Vertex_handle(3)] = 9;
+	expected_num_simplices[Vertex_handle(4)] = 7;
+	expected_num_simplices[Vertex_handle(5)] = 7;
+
+	for(auto pair : expected_num_simplices){
+		unsigned num_simplices_around = 0;
+		for(const auto& simplex : complex.simplex_range(pair.first)){
+			simplex.dimension();
+			DBGVALUE(simplex);
+			++num_simplices_around;
+		}
+
+		correct_number_simplices = correct_number_simplices && (num_simplices_around == pair.second);
+
+		DBGMSG("current vertex:",pair.first);
+		DBGMSG("expected_num_simplices:",pair.second);
+		DBGMSG("found:",num_simplices_around);
+	}
+	return correct_number_simplices;
+}
+
+
+
+bool test_iterator_simplices2(){
+	Complex complex(2);
+	complex.add_edge(Vertex_handle(0),Vertex_handle(1));
+
+	for(const auto& s:complex.triangle_range()){
+		s.dimension();
+		return false; // there are no triangles
+	}
+
+	unsigned num_simplices = 0 ;
+
+
+	DBGVALUE(complex.to_string());
+
+	for(const auto& simplex : complex.simplex_range(Vertex_handle(0))){
+		simplex.dimension();
+		DBGVALUE(simplex);
+	}
+
+
+	for(const auto& simplex : complex.simplex_range()){
+		DBGVALUE(simplex);
+		simplex.dimension();
+		++num_simplices;
+	}
+	bool correct_number_simplices = (num_simplices == 3);
+	return correct_number_simplices;
+}
+
+
+bool test_iterator_simplices3(){
+	Complex complex(3);
+	complex.add_edge(Vertex_handle(0),Vertex_handle(1));
+	complex.add_edge(Vertex_handle(1),Vertex_handle(2));
+	complex.add_edge(Vertex_handle(2),Vertex_handle(0));
+	complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2));
+
+	unsigned num_simplices = 0 ;
+
+	for(const auto& simplex : complex.simplex_range(Vertex_handle(0))){
+		simplex.dimension();
+		DBGVALUE(simplex);
+	}
+
+
+	for(const auto& simplex : complex.simplex_range()){
+		DBGVALUE(simplex);
+		simplex.dimension();
+		++num_simplices;
+	}
+	bool correct_number_simplices = (num_simplices == 6);
+	return correct_number_simplices;
+}
+
+bool test_iterator_simplices4(){
+	Complex empty_complex;
+	for(auto v : empty_complex.vertex_range()){
+		v;
+	}
+	for(auto e : empty_complex.edge_range()){
+		empty_complex[e];
+	}
+	for(auto t : empty_complex.triangle_range()){
+		t.dimension();
+	}
+	for(auto s : empty_complex.simplex_range()){
+		s.dimension();
+	}
+	return true;
+}
+
+
+
+
+
+template<typename Map>
+auto blocker_range(Map map) -> decltype( map | boost::adaptors::map_values){
+	return map| boost::adaptors::map_values ;
+}
+
+
+bool test_iterator_blockers(){
+	Complex complex;
+	Simplex_handle alpha;
+	Simplex_handle vertex_set_expected;
+	// Build the complexes
+	for (int i=0;i<20;i++){
+		complex.add_vertex();
+	}
+	for (int i=10;i<15;i++){
+		for (int j=i+1;j<15;j++)
+			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+	}
+
+	complex.add_blocker(Simplex_handle(Vertex_handle(10),Vertex_handle(11),Vertex_handle(12)));
+	complex.add_blocker(Simplex_handle(Vertex_handle(2),Vertex_handle(1),Vertex_handle(10)));
+	complex.add_blocker(Simplex_handle(Vertex_handle(10),Vertex_handle(9),Vertex_handle(15)));
+	complex.add_blocker(Simplex_handle(Vertex_handle(1),Vertex_handle(9),Vertex_handle(8)));
+
+	// Print result
+	int num_blockers=0;
+	for(auto blockers : complex.blocker_range(Vertex_handle(10))){
+		TESTVALUE(*blockers) ;
+		num_blockers++;
+	}
+	bool test = (num_blockers==3);
+
+	num_blockers=0;
+	for (auto blockers : complex.blocker_range()){
+		TESTVALUE(*blockers) ;
+		num_blockers++;
+	}
+	test = test && (num_blockers==4) ;
+
+	return test;
+}
+
+
+bool test_link0(){
+
+	enum { a, b, c, d, n };
+	Complex complex(n);
+	complex.add_edge(Vertex_handle(b),Vertex_handle(c));complex.add_edge(Vertex_handle(c),Vertex_handle(d));
+	Simplex_handle alpha = Simplex_handle(Vertex_handle(c));
+	Skeleton_blocker_link_complex<Complex> L(complex,alpha);
+	PRINT(L.num_vertices());
+	PRINT(L.to_string());
+
+	bool test1 = L.contains_vertex(*L.get_address(Root_vertex_handle(b)));
+	bool test2 = L.contains_vertex(*L.get_address(Root_vertex_handle(d)));
+	bool test3 = L.num_edges()==0;
+	bool test4 = L.num_blockers()==0;
+	return test1&&test2&&test3&&test4;
+
+}
+
+bool test_link1(){
+	Complex complex;
+
+
+	// Build the complexes
+	for (int i=0;i<20;i++){
+		complex.add_vertex();
+	}
+	for (int i=10;i<15;i++){
+		for (int j=i+1;j<15;j++)
+			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+	}
+	Simplex_handle alpha(Vertex_handle(12),Vertex_handle(14));
+	Skeleton_blocker_link_complex<Complex> L(complex,alpha);
+	// Complexes built
+
+	// verification
+	bool test1 = L.contains_vertex(*L.get_address(Root_vertex_handle(10)));
+	bool test2 = L.contains_vertex(*L.get_address(Root_vertex_handle(11)));
+	bool test3 = L.contains_vertex(*L.get_address(Root_vertex_handle(13)));
+	bool test4 = L.num_edges()==3;
+	bool test5 = L.num_blockers()==0;
+	Root_simplex_handle simplex;
+	simplex.add_vertex(Root_vertex_handle(10));
+	simplex.add_vertex(Root_vertex_handle(11));
+	simplex.add_vertex(Root_vertex_handle(13));
+	bool test6(L.get_simplex_address(simplex));
+	bool test7 = L.contains(*(L.get_simplex_address(simplex)));
+	cerr <<"----> Ocomplex \n";
+	return test1&&test2&&test3&&test4&&test5&&test6&&test7 ;
+
+}
+
+
+bool test_link2(){
+	Complex complex;
+
+	Simplex_handle alpha;
+	Simplex_handle vertex_set_expected;
+	// Build the complexes
+	for (int i=0;i<20;i++){
+		complex.add_vertex();
+	}
+	for (int i=10;i<15;i++){
+		for (int j=i+1;j<15;j++)
+			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+	}
+	complex.add_blocker(Vertex_handle(10),Vertex_handle(11),Vertex_handle(13));
+	alpha = Simplex_handle(Vertex_handle(12),Vertex_handle(14));
+	Skeleton_blocker_link_complex<Complex> L(complex,alpha);
+	// Complexes built
+
+	// Print result
+	cerr << "complex complex"<< complex.to_string();
+	cerr <<endl<<endl;
+	cerr << "L= Link_complex("<<alpha<<") : \n"<<L.to_string();
+
+	// verification
+	bool test1 = L.contains_vertex(*L.get_address(Root_vertex_handle(10)));
+	bool test2 = L.contains_vertex(*L.get_address(Root_vertex_handle(11)));
+	bool test3 = L.contains_vertex(*L.get_address(Root_vertex_handle(13)));
+	bool test4 = L.num_edges()==3;
+	bool test5 = L.num_blockers()==1;
+	Root_simplex_handle simplex;
+	simplex.add_vertex(Root_vertex_handle(10));
+	simplex.add_vertex(Root_vertex_handle(11));
+	simplex.add_vertex(Root_vertex_handle(13));
+	bool test6 = L.contains_blocker(*(L.get_simplex_address(simplex)));
+	cerr <<"----> Ocomplex \n";
+	return test1&&test2&&test3&&test4&&test5&&test6 ;
+}
+
+bool test_link3(){
+	Complex complex;
+
+	Simplex_handle alpha;
+	Simplex_handle vertex_set_expected;
+	// Build the complexes
+	for (int i=0;i<20;i++){
+		complex.add_vertex();
+	}
+	for (int i=10;i<15;i++){
+		for (int j=i+1;j<15;j++)
+			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+	}
+	complex.add_blocker(Vertex_handle(10),Vertex_handle(11),Vertex_handle(12));
+	alpha = Simplex_handle(Vertex_handle(12),Vertex_handle(14));
+	Skeleton_blocker_link_complex<Complex> L(complex,alpha);
+	// Complexes built
+
+	// Print result
+	cerr << "complex complex"<< complex.to_string();
+	cerr <<endl<<endl;
+	cerr << "L= Link_complex("<<alpha<<") : \n"<<L.to_string();
+
+	// verification
+	bool test = assert_vertex(L,*L.get_address(Root_vertex_handle(10)));
+	test = test&& assert_vertex(L,*L.get_address(Root_vertex_handle(11)));
+	test = test&& assert_vertex(L,*L.get_address(Root_vertex_handle(13)));
+	test = test&& L.num_edges()==2;
+	test = test&&L.contains_edge(*L.get_address(Root_vertex_handle(10)),*L.get_address(Root_vertex_handle(13)));
+	test=test&&L.contains_edge(*L.get_address(Root_vertex_handle(13)),*L.get_address(Root_vertex_handle(11)));
+	test=test&&L.num_blockers()==0;
+	return test;
+}
+
+bool test_link4(){
+	Complex complex;
+
+	// Build the complexes
+	for (int i=0;i<20;i++){
+		complex.add_vertex();
+	}
+	for (int i=10;i<15;i++){
+		for (int j=i+1;j<15;j++)
+			complex.add_edge(Vertex_handle(i),Vertex_handle(j));
+	}
+	complex.add_blocker(Vertex_handle(10),Vertex_handle(11),Vertex_handle(12),Vertex_handle(13));
+	Simplex_handle alpha(Vertex_handle(12),Vertex_handle(14));
+	Skeleton_blocker_link_complex<Complex> L(complex,alpha);
+	// Complexes built
+
+	// verification
+	bool test1 = L.contains_vertex(*L.get_address(Root_vertex_handle(10)));
+	bool test2 = L.contains_vertex(*L.get_address(Root_vertex_handle(11)));
+	bool test3 = L.contains_vertex(*L.get_address(Root_vertex_handle(13)));
+	bool test4 = L.num_edges()==3;
+	bool test5 = L.num_blockers()==1;
+	Root_simplex_handle simplex;
+	simplex.add_vertex(Root_vertex_handle(10));
+	simplex.add_vertex(Root_vertex_handle(11));
+	simplex.add_vertex(Root_vertex_handle(13));
+	bool test6 = L.contains_blocker(*(L.get_simplex_address(simplex)));
+	cerr <<"----> Ocomplex \n";
+	return test1&&test2&&test3&&test4&&test5&&test6 ;
+
+}
+
+bool test_link5(){
+	Complex complex(0,new Print_complex_visitor<Vertex_handle>());
+	// Build the complexes
+	build_complete(4,complex);
+	complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2),Vertex_handle(3));
+
+	Simplex_handle alpha(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2));
+
+
+	Skeleton_blocker_link_complex<Complex> L(complex,alpha);	// Complexes built
+
+	// Print result
+	PRINT(complex.to_string());
+	cerr <<endl<<endl;
+	PRINT(L.to_string());
+
+	// verification
+	return L.num_vertices()==0;
+}
+
+bool test_link6(){
+	Complex complex(0,new Print_complex_visitor<Vertex_handle>());
+	// Build the complexes
+	build_complete(4,complex);
+	complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2));
+
+	Simplex_handle alpha(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2));
+
+	Skeleton_blocker_link_complex<Complex> link_blocker_alpha;
+
+	build_link_of_blocker(complex,alpha,link_blocker_alpha);
+
+	// Print result
+	PRINT(complex.to_string());
+	cerr <<endl<<endl;
+	PRINT(link_blocker_alpha.to_string());
+
+	// verification
+	return link_blocker_alpha.num_vertices()==1;
+}
+
+
+bool test_link7(){
+	Complex complex(0,new Print_complex_visitor<Vertex_handle>());
+	// Build the complexes
+	build_complete(6,complex);
+	complex.add_vertex();
+	complex.add_vertex();
+	for(int i = 3; i<6; ++i){
+		complex.add_edge(Vertex_handle(i),Vertex_handle(6));
+		complex.add_edge(Vertex_handle(i),Vertex_handle(7));
+	}
+	complex.add_edge(Vertex_handle(6),Vertex_handle(7));
+	complex.add_blocker(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2));
+	complex.add_blocker(Vertex_handle(3),Vertex_handle(4),Vertex_handle(5));
+
+	Simplex_handle alpha(Vertex_handle(3),Vertex_handle(4),Vertex_handle(5));
+
+	Skeleton_blocker_link_complex<Complex> link_blocker_alpha;
+
+	build_link_of_blocker(complex,alpha,link_blocker_alpha);
+
+	//the result should be the edge {6,7} plus the blocker {0,1,2}
+
+	// Print result
+	PRINT(complex.to_string());
+	cerr <<endl<<endl;
+	DBGVALUE(link_blocker_alpha.to_string());
+
+	Skeleton_blocker_link_complex<Complex> link_blocker_alpha_cpy = link_blocker_alpha;
+
+	DBGVALUE(link_blocker_alpha_cpy.to_string());
+
+	bool equal_complexes =
+			(link_blocker_alpha.num_vertices() == link_blocker_alpha_cpy.num_vertices())
+			&&(link_blocker_alpha.num_blockers() == link_blocker_alpha_cpy.num_blockers())
+			&&(link_blocker_alpha.num_edges() == link_blocker_alpha_cpy.num_edges())
+			;
+	DBGVALUE((link_blocker_alpha.num_blockers() == link_blocker_alpha_cpy.num_blockers()));
+	DBGVALUE((link_blocker_alpha.num_blockers() ));
+	DBGVALUE(( link_blocker_alpha_cpy.num_blockers()));
+
+	DBGVALUE(equal_complexes);
+
+	// verification
+	return link_blocker_alpha.num_vertices()==5 && link_blocker_alpha.num_edges()==4 && link_blocker_alpha.num_blockers()==1 && equal_complexes;
+}
+
+
+
+
+
+
+
+
+template<typename SimplexHandle>
+void add_triangle_edges(int a,int b,int c,list<SimplexHandle>& simplices){
+	typedef SimplexHandle Simplex_handle;
+	typedef typename SimplexHandle::Vertex_handle Vertex_handle;
+
+	simplices.push_back(Simplex_handle(Vertex_handle(a),Vertex_handle(b) ));
+	simplices.push_back(Simplex_handle(Vertex_handle(b),Vertex_handle(c) ));
+	simplices.push_back(Simplex_handle(Vertex_handle(c),Vertex_handle(a) ));
+}
+
+template<typename SimplexHandle>
+void add_triangle(int a,int b,int c,list<SimplexHandle>& simplices){
+	typedef SimplexHandle Simplex_handle;
+	typedef typename SimplexHandle::Vertex_handle Vertex_handle;
+	simplices.push_back(Simplex_handle(Vertex_handle(a),Vertex_handle(b),Vertex_handle(c)));
+}
+
+bool test_constructor(){
+	list <Simplex_handle> simplices;
+
+	simplices.push_back(Simplex_handle(Vertex_handle(0)));
+	simplices.push_back(Simplex_handle(Vertex_handle(1)));
+	simplices.push_back(Simplex_handle(Vertex_handle(2)));
+	simplices.push_back(Simplex_handle(Vertex_handle(3)));
+	simplices.push_back(Simplex_handle(Vertex_handle(4)));
+	simplices.push_back(Simplex_handle(Vertex_handle(5)));
+
+	simplices.push_back(Simplex_handle(Vertex_handle(3),Vertex_handle(5) ));
+
+	add_triangle_edges(0,1,5,simplices);
+	add_triangle_edges(1,2,3,simplices);
+	add_triangle_edges(2,3,4,simplices);
+	add_triangle_edges(1,3,4,simplices);
+	add_triangle_edges(1,2,4,simplices);
+
+
+	add_triangle(0,1,5,simplices);
+	add_triangle(1,2,3,simplices);
+	add_triangle(2,3,4,simplices);
+	add_triangle(1,3,4,simplices);
+	add_triangle(1,2,4,simplices);
+
+
+	Complex complex(simplices);
+
+	PRINT(complex.to_string());
+
+	return ( complex.num_vertices()==6&&complex.num_edges()==10&&  complex.num_blockers()==2);
+}
+
+
+list<Simplex_handle> subfaces(Simplex_handle top_face){
+	list<Simplex_handle> res;
+	if(top_face.dimension()==-1) return res;
+	if(top_face.dimension()==0) {
+		res.push_back(top_face);
+		return res;
+	}
+	else{
+		Vertex_handle first_vertex = top_face.first_vertex();
+		top_face.remove_vertex(first_vertex);
+		res = subfaces(top_face);
+		list<Simplex_handle> copy = res;
+		for(auto& simplex : copy){
+			simplex.add_vertex(first_vertex);
+		}
+		res.push_back(Simplex_handle(first_vertex));
+		res.splice(res.end(),copy);
+		return res;
+	}
+}
+
+
+bool test_constructor2(){
+	Simplex_handle simplex;
+	for(int i =0 ; i < 5;++i)
+		simplex.add_vertex(i);
+
+	list <Simplex_handle> simplices(subfaces(simplex));
+	simplices.remove(simplex);
+
+	Complex complex(simplices);
+
+	PRINT(complex.to_string());
+
+	for(auto b : complex.const_blocker_range()){
+		cout << "b:"<<b<<endl;
+	}
+
+	return ( complex.num_vertices()==5&&complex.num_edges()==10&&  complex.num_blockers()==1);
+}
+
+
+
+
+int main (int argc, char *argv[])
+{
+	Tests tests_complex;
+	tests_complex.add("test simplex",test_simplex);
+	tests_complex.add("test_link0",test_link0);
+	tests_complex.add("test_link1",test_link1);
+	tests_complex.add("test_link2",test_link2);
+	tests_complex.add("test_link3",test_link3);
+	tests_complex.add("test_link4",test_link4);
+	tests_complex.add("test_link5",test_link5);
+	tests_complex.add("test_link6",test_link6);
+	tests_complex.add("test_link7",test_link7);
+
+	tests_complex.add("test iterator vertices 1",test_iterator_vertices1);
+	tests_complex.add("test iterator vertices 2",test_iterator_vertices2);
+	tests_complex.add("test iterator edges",test_iterator_edge);
+	tests_complex.add("test iterator edges 2",test_iterator_edge2);
+	tests_complex.add("test iterator edges 3",test_iterator_edge3);
+
+	tests_complex.add("test iterator simplices",test_iterator_simplices);
+	tests_complex.add("test iterator simplices2",test_iterator_simplices2);
+	tests_complex.add("test iterator simplices3",test_iterator_simplices3);
+	tests_complex.add("test iterator simplices4",test_iterator_simplices4);
+
+
+	tests_complex.add("test iterator blockers",test_iterator_blockers);
+	tests_complex.add("test_iterator_triangles",test_iterator_triangles);
+
+	tests_complex.add("test_constructor_list_simplices",test_constructor);
+	tests_complex.add("test_constructor_list_simplices2",test_constructor2);
+
+	if(tests_complex.run()){
+		return EXIT_SUCCESS;
+	}
+	else{
+		return EXIT_FAILURE;
+	}
+
+	//	test_iterator_simplices();
+}
+
-- 
cgit v1.2.3