Merged latest trunk to witness

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

Former-commit-id: a90ff18c8991da76cf66daec3dc8f68c049cfe57

32 files changed, 4078 insertions, 2593 deletions

diff --git a/src/Alpha_shapes/example/CMakeLists.txt b/src/Alpha_shapes/example/CMakeLists.txt
new file mode 100644
index 00000000..fb94ca05
--- /dev/null
+++ b/src/Alpha_shapes/example/CMakeLists.txt
@@ -0,0 +1,31 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIAlphaShapesExample)
+
+# need CGAL 4.6
+# cmake -DCGAL_DIR=~/workspace/CGAL-4.6-beta1 ../../..
+if(CGAL_FOUND)
+  if (NOT CGAL_VERSION VERSION_LESS 4.6.0)
+    message(STATUS "CGAL version: ${CGAL_VERSION}.")
+
+    include( ${CGAL_USE_FILE} )
+ 
+    find_package(Eigen3 3.1.0)
+    if (EIGEN3_FOUND)
+      message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.")
+      include( ${EIGEN3_USE_FILE} )
+      include_directories (BEFORE "../../include")
+
+      add_executable ( dtoffrw Delaunay_triangulation_off_rw.cpp )
+      target_link_libraries(dtoffrw  ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY})
+      add_test(dtoffrw_tore3D ${CMAKE_CURRENT_BINARY_DIR}/dtoffrw ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off 3)
+
+      #add_definitions(-DDEBUG_TRACES)
+      add_executable ( stfromdt Simplex_tree_from_delaunay_triangulation.cpp )
+      target_link_libraries(stfromdt  ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY})
+    else()
+      message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.")
+    endif()
+  else()
+    message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.")
+  endif ()
+endif()
diff --git a/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp
new file mode 100644
index 00000000..03f6440d
--- /dev/null
+++ b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp
@@ -0,0 +1,105 @@
+/*    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):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2014  INRIA Saclay (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/>.
+ */
+
+// to construct a Delaunay_triangulation from a OFF file
+#include "gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h"
+
+#include <CGAL/Delaunay_triangulation.h>
+#include <CGAL/Epick_d.h>
+#include <CGAL/point_generators_d.h>
+#include <CGAL/algorithm.h>
+#include <CGAL/assertions.h>
+
+#include <iostream>
+#include <iterator>
+#include <vector>
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+
+// Use dynamic_dimension_tag for the user to be able to set dimension
+typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > K;
+typedef CGAL::Delaunay_triangulation<K> T;
+// The triangulation uses the default instanciation of the 
+// TriangulationDataStructure template parameter
+
+void usage(char * const progName) {
+  std::cerr << "Usage: " << progName << " filename.off dimension" << std::endl;
+  exit(-1); // ----- >>
+}
+
+int main(int argc, char **argv) {
+  if (argc != 3) {
+    std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl;
+    usage(argv[0]);
+  }
+
+  int dimension = 0;
+  int returnedScanValue = sscanf(argv[2], "%d", &dimension);
+  if ((returnedScanValue == EOF) || (dimension <= 0)) {
+    std::cerr << "Error: " << argv[2] << " is not correct" << std::endl;
+    usage(argv[0]);
+  }
+
+  T dt(dimension);
+  std::string offFileName(argv[1]);
+  Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(offFileName, dt, true, true);
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read file " << offFileName << std::endl;
+    exit(-1); // ----- >>
+  }
+
+  std::cout << "number of vertices=" << dt.number_of_vertices() << std::endl;
+  std::cout << "number of full cells=" << dt.number_of_full_cells() << std::endl;
+  std::cout << "number of finite full cells=" << dt.number_of_finite_full_cells() << std::endl;
+
+  // Points list
+  /*for (T::Vertex_iterator vit = dt.vertices_begin(); vit != dt.vertices_end(); ++vit) {
+        for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) {
+          std::cout << *Coord << " ";
+        }
+        std::cout << std::endl;
+  }
+  std::cout << std::endl;*/
+ 
+  int i = 0, j = 0;
+  typedef T::Full_cell_iterator Full_cell_iterator;
+  typedef T::Facet Facet;
+
+  for (Full_cell_iterator cit = dt.full_cells_begin(); cit != dt.full_cells_end(); ++cit) {
+    if (!dt.is_infinite(cit)) {
+      j++;
+      continue;
+    }
+    Facet fct(cit, cit->index(dt.infinite_vertex()));
+    i++;
+  }
+  std::cout << "There are " << i << " facets on the convex hull." << std::endl;
+  std::cout << "There are " << j << " facets not on the convex hull." << std::endl;
+
+
+  std::string offOutputFile("out.off");
+  Gudhi::alphashapes::Delaunay_triangulation_off_writer<T> off_writer(offOutputFile, dt);
+
+  return 0;
+}
\ No newline at end of file
diff --git a/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp b/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp
new file mode 100644
index 00000000..3a17b519
--- /dev/null
+++ b/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp
@@ -0,0 +1,103 @@
+/*    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):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2014  INRIA Saclay (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/>.
+ */
+
+// to construct a Delaunay_triangulation from a OFF file
+#include "gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h"
+#include "gudhi/Alpha_shapes.h"
+
+// to construct a simplex_tree from Delaunay_triangulation
+#include "gudhi/graph_simplicial_complex.h"
+#include "gudhi/Simplex_tree.h"
+
+#include <CGAL/Delaunay_triangulation.h>
+#include <CGAL/Epick_d.h>
+#include <CGAL/point_generators_d.h>
+#include <CGAL/algorithm.h>
+#include <CGAL/assertions.h>
+
+#include <iostream>
+#include <iterator>
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+
+// Use dynamic_dimension_tag for the user to be able to set dimension
+typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > K;
+typedef CGAL::Delaunay_triangulation<K> T;
+// The triangulation uses the default instanciation of the 
+// TriangulationDataStructure template parameter
+
+void usage(char * const progName) {
+  std::cerr << "Usage: " << progName << " filename.off dimension" << std::endl;
+  exit(-1); // ----- >>
+}
+
+int main(int argc, char **argv) {
+  if (argc != 3) {
+    std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl;
+    usage(argv[0]);
+  }
+
+  int dimension = 0;
+  int returnedScanValue = sscanf(argv[2], "%d", &dimension);
+  if ((returnedScanValue == EOF) || (dimension <= 0)) {
+    std::cerr << "Error: " << argv[2] << " is not correct" << std::endl;
+    usage(argv[0]);
+  }
+
+  // ----------------------------------------------------------------------------
+  //
+  // Init of an alpha-shape from a Delaunay triangulation
+  //
+  // ----------------------------------------------------------------------------
+  T dt(dimension);
+  std::string off_file_name(argv[1]);
+
+  Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(off_file_name, dt, false, false);
+  if (!off_reader.is_valid()) {
+    std::cerr << "Unable to read file " << off_file_name << std::endl;
+    exit(-1); // ----- >>
+  }
+
+  std::cout << "number of vertices=" << dt.number_of_vertices() << std::endl;
+  std::cout << "number of full cells=" << dt.number_of_full_cells() << std::endl;
+  std::cout << "number of finite full cells=" << dt.number_of_finite_full_cells() << std::endl;
+
+  Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_dt(dt);
+  //std::cout << alpha_shapes_from_dt << std::endl;
+
+  // ----------------------------------------------------------------------------
+  //
+  // Init of an alpha-shape from a OFF file
+  //
+  // ----------------------------------------------------------------------------
+  Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_file(off_file_name, dimension);
+  //std::cout << alpha_shapes_from_file << std::endl;
+  
+  std::cout << "alpha_shapes_from_file.dimension()=" << alpha_shapes_from_file.dimension() << std::endl;
+  std::cout << "alpha_shapes_from_file.filtration()=" << alpha_shapes_from_file.filtration() << std::endl;
+  std::cout << "alpha_shapes_from_file.num_simplices()=" << alpha_shapes_from_file.num_simplices() << std::endl;
+  std::cout << "alpha_shapes_from_file.num_vertices()=" << alpha_shapes_from_file.num_vertices() << std::endl;
+
+  return 0;
+}
\ No newline at end of file
diff --git a/src/Alpha_shapes/include/gudhi/Alpha_shapes.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h
new file mode 100644
index 00000000..b8efdb4d
--- /dev/null
+++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h
@@ -0,0 +1,200 @@
+/*    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):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_
+#define SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_
+
+// to construct a Delaunay_triangulation from a OFF file
+#include <gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h>
+
+// to construct a simplex_tree from Delaunay_triangulation
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/Simplex_tree.h>
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include <CGAL/Delaunay_triangulation.h>
+#include <CGAL/Epick_d.h>
+#include <CGAL/algorithm.h>
+#include <CGAL/assertions.h>
+
+#include <iostream>
+#include <iterator>
+#include <vector>
+#include <string>
+
+namespace Gudhi {
+
+namespace alphashapes {
+
+/** \defgroup alpha_shapes Alpha shapes in dimension N
+ *
+ <DT>Implementations:</DT>
+ Alpha shapes in dimension N are a subset of Delaunay Triangulation in dimension N.
+
+
+ * \author    Vincent Rouvreau
+ * \version   1.0
+ * \date      2015
+ * \copyright GNU General Public License v3.
+ * @{
+ */
+
+/**
+ * \brief Alpha shapes data structure.
+ *
+ * \details Every simplex \f$[v_0, \cdots ,v_d]\f$ admits a canonical orientation
+ * induced by the order relation on vertices \f$ v_0 < \cdots < v_d \f$.
+ *
+ * Details may be found in \cite boissonnatmariasimplextreealgorithmica.
+ *
+ * \implements FilteredComplex
+ *
+ */
+class Alpha_shapes {
+ private:
+  // From Simplex_tree
+  /** \brief Type required to insert into a simplex_tree (with or without subfaces).*/
+  typedef std::vector<Vertex_handle> typeVectorVertex;
+
+  // From CGAL
+  /** \brief Kernel for the Delaunay_triangulation.
+   * Dimension can be set dynamically.
+   */
+  typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel;
+  /** \brief Delaunay_triangulation type required to create an alpha-shape.
+   */
+  typedef CGAL::Delaunay_triangulation<Kernel> Delaunay_triangulation;
+
+ private:
+  /** \brief Upper bound on the simplex tree of the simplicial complex.*/
+  Gudhi::Simplex_tree<> _st;
+
+ public:
+
+  Alpha_shapes(std::string off_file_name, int dimension) {
+    Delaunay_triangulation dt(dimension);
+    Gudhi::alphashapes::Delaunay_triangulation_off_reader<Delaunay_triangulation>
+        off_reader(off_file_name, dt, true, true);
+    if (!off_reader.is_valid()) {
+      std::cerr << "Unable to read file " << off_file_name << std::endl;
+      exit(-1); // ----- >>
+    }
+#ifdef DEBUG_TRACES
+    std::cout << "number of vertices=" << dt.number_of_vertices() << std::endl;
+    std::cout << "number of full cells=" << dt.number_of_full_cells() << std::endl;
+    std::cout << "number of finite full cells=" << dt.number_of_finite_full_cells() << std::endl;
+#endif  // DEBUG_TRACES
+    init<Delaunay_triangulation>(dt);
+  }
+
+  template<typename T>
+  Alpha_shapes(T triangulation) {
+    init<T>(triangulation);
+  }
+
+  ~Alpha_shapes() { }
+
+ private:
+
+  template<typename T>
+  void init(T triangulation) {
+    _st.set_dimension(triangulation.maximal_dimension());
+    _st.set_filtration(0.0);
+    // triangulation points list
+    for (auto vit = triangulation.finite_vertices_begin();
+         vit != triangulation.finite_vertices_end(); ++vit) {
+      typeVectorVertex vertexVector;
+      Vertex_handle vertexHdl = std::distance(triangulation.finite_vertices_begin(), vit);
+      vertexVector.push_back(vertexHdl);
+
+      // Insert each point in the simplex tree
+      _st.insert_simplex(vertexVector, 0.0);
+
+#ifdef DEBUG_TRACES
+      std::cout << "P" << vertexHdl << ":";
+      for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) {
+        std::cout << *Coord << " ";
+      }
+      std::cout << std::endl;
+#endif  // DEBUG_TRACES
+    }
+    // triangulation finite full cells list
+    for (auto cit = triangulation.finite_full_cells_begin();
+         cit != triangulation.finite_full_cells_end(); ++cit) {
+      typeVectorVertex vertexVector;
+      for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
+        // Vertex handle is distance - 1
+        Vertex_handle vertexHdl = std::distance(triangulation.vertices_begin(), *vit) - 1;
+        vertexVector.push_back(vertexHdl);
+      }
+      // Insert each point in the simplex tree
+      _st.insert_simplex_and_subfaces(vertexVector, 0.0);
+
+#ifdef DEBUG_TRACES
+      std::cout << "C" << std::distance(triangulation.finite_full_cells_begin(), cit) << ":";
+      for (auto value : vertexVector) {
+        std::cout << value << ' ';
+      }
+      std::cout << std::endl;
+#endif  // DEBUG_TRACES
+    }
+  }
+
+ public:
+
+  /** \brief Returns the number of vertices in the complex. */
+  size_t num_vertices() {
+    return _st.num_vertices();
+  }
+
+  /** \brief Returns the number of simplices in the complex.
+   *
+   * Does not count the empty simplex. */
+  const unsigned int& num_simplices() const {
+    return _st.num_simplices();
+  }
+
+  /** \brief Returns an upper bound on the dimension of the simplicial complex. */
+  int dimension() {
+    return _st.dimension();
+  }
+
+  /** \brief Returns an upper bound of the filtration values of the simplices. */
+  Filtration_value filtration() {
+    return _st.filtration();
+  }
+
+  friend std::ostream& operator<<(std::ostream& os, const Alpha_shapes& alpha_shape) {
+    // TODO: Program terminated with signal SIGABRT, Aborted - Maybe because of copy constructor
+    Gudhi::Simplex_tree<> st = alpha_shape._st;
+    os << st << std::endl;
+    return os;
+  }
+};
+
+} // namespace alphashapes
+
+} // namespace Gudhi
+
+#endif  // SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_H_
diff --git a/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h
new file mode 100644
index 00000000..693b393e
--- /dev/null
+++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h
@@ -0,0 +1,213 @@
+/*    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):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#ifndef SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_
+#define SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_
+
+#include <string>
+#include <vector>
+#include <fstream>
+
+#include "gudhi/Off_reader.h"
+
+namespace Gudhi {
+
+namespace alphashapes {
+
+/**
+ *@brief Off reader visitor with flag that can be passed to Off_reader to read a Delaunay_triangulation_complex.
+ */
+template<typename Complex>
+class Delaunay_triangulation_off_flag_visitor_reader {
+  Complex& complex_;
+  typedef typename Complex::Point Point;
+
+  const bool load_only_points_;
+
+ public:
+  explicit Delaunay_triangulation_off_flag_visitor_reader(Complex& complex, bool load_only_points = false) :
+      complex_(complex),
+      load_only_points_(load_only_points) { }
+
+  void init(int dim, int num_vertices, int num_faces, int num_edges) {
+#ifdef DEBUG_TRACES
+    std::cout << "init" << std::endl;
+#endif  // DEBUG_TRACES
+  }
+
+  void point(const std::vector<double>& point) {
+#ifdef DEBUG_TRACES
+    std::cout << "p ";
+    for (auto coordinate: point) {
+      std::cout << coordinate << " | ";
+    }
+    std::cout << std::endl;
+#endif  // DEBUG_TRACES
+    complex_.insert(Point(point.size(), point.begin(), point.end()));
+  }
+
+  void maximal_face(const std::vector<int>& face) {
+    // For alpha shapes, only points are read
+  }
+
+  void done() {
+#ifdef DEBUG_TRACES
+    std::cout << "done" << std::endl;
+#endif  // DEBUG_TRACES
+  }
+};
+
+/**
+ *@brief Off reader visitor that can be passed to Off_reader to read a Delaunay_triangulation_complex.
+ */
+template<typename Complex>
+class Delaunay_triangulation_off_visitor_reader {
+  Complex& complex_;
+  // typedef typename Complex::Vertex_handle Vertex_handle;
+  // typedef typename Complex::Simplex_handle Simplex_handle;
+  typedef typename Complex::Point Point;
+
+  const bool load_only_points_;
+  std::vector<Point> points_;
+  // std::vector<Simplex_handle> maximal_faces_;
+
+ public:
+  explicit Delaunay_triangulation_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) {
+#ifdef DEBUG_TRACES
+    std::cout << "init - " << num_vertices << std::endl;
+#endif  // DEBUG_TRACES
+    // maximal_faces_.reserve(num_faces);
+    points_.reserve(num_vertices);
+  }
+
+  void point(const std::vector<double>& point) {
+#ifdef DEBUG_TRACES
+    std::cout << "p ";
+    for (auto coordinate: point) {
+      std::cout << coordinate << " | ";
+    }
+    std::cout << std::endl;
+#endif  // DEBUG_TRACES
+    points_.emplace_back(Point(point.size(), point.begin(), point.end()));
+  }
+
+  void maximal_face(const std::vector<int>& face) {
+    // For alpha shapes, only points are read
+  }
+
+  void done() {
+    complex_.insert(points_.begin(), points_.end());
+#ifdef DEBUG_TRACES
+    std::cout << "done" << std::endl;
+#endif  // DEBUG_TRACES
+  }
+};
+
+/**
+ *@brief Class that allows to load a Delaunay_triangulation_complex from an off file.
+ */
+template<typename Complex>
+class Delaunay_triangulation_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
+   */
+  Delaunay_triangulation_off_reader(const std::string & name_file, Complex& read_complex, bool read_only_points = false,
+                                    bool is_flag = false) : valid_(false) {
+    std::ifstream stream(name_file);
+    if (stream.is_open()) {
+      if (is_flag) {
+        // For alpha shapes, only points are read
+        Delaunay_triangulation_off_flag_visitor_reader<Complex> off_visitor(read_complex, true);
+        Off_reader off_reader(stream);
+        valid_ = off_reader.read(off_visitor);
+      } else {
+        // For alpha shapes, only points are read
+        Delaunay_triangulation_off_visitor_reader<Complex> off_visitor(read_complex, true);
+        Off_reader off_reader(stream);
+        valid_ = off_reader.read(off_visitor);
+      }
+    }
+  }
+
+  /**
+   * return true if reading did not meet problems.
+   */
+  bool is_valid() const {
+    return valid_;
+  }
+
+ private:
+  bool valid_;
+};
+
+template<typename Complex>
+class Delaunay_triangulation_off_writer {
+ public:
+  /**
+   * name_file : file where the off will be written
+   * save_complex : complex that be outputted in the file
+   * for now only save triangles.
+   */
+  Delaunay_triangulation_off_writer(const std::string & name_file, const Complex& save_complex) {
+    std::ofstream stream(name_file);
+    if (stream.is_open()) {
+      // OFF header
+      stream << "OFF" << std::endl;
+      // no endl on next line - don't know why...
+      stream << save_complex.number_of_vertices() << " " << save_complex.number_of_finite_full_cells() << " 0";
+
+      // Points list
+      for (auto vit = save_complex.vertices_begin(); vit != save_complex.vertices_end(); ++vit) {
+        for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) {
+          stream << *Coord << " ";
+        }
+        stream << std::endl;
+      }
+
+      // Finite cells list
+      for (auto cit = save_complex.finite_full_cells_begin(); cit != save_complex.finite_full_cells_end(); ++cit) {
+        stream << std::distance(cit->vertices_begin(), cit->vertices_end()) << " ";  // Dimension
+        for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
+          auto vertexHdl = *vit;
+          // auto vertexHdl = std::distance(save_complex.vertices_begin(), *vit) - 1;
+          // stream << std::distance(save_complex.vertices_begin(), *(vit)) - 1 << " ";
+        }
+        stream << std::endl;
+      }
+      stream.close();
+    } else {
+      std::cerr << "could not open file " << name_file << std::endl;
+    }
+  }
+};
+
+}  // namespace alphashapes
+
+}  // namespace Gudhi
+
+#endif  // SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_
diff --git a/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp b/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp
new file mode 100644
index 00000000..a90704b6
--- /dev/null
+++ b/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp
@@ -0,0 +1,126 @@
+/*    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):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2015  INRIA Saclay (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/>.
+ */
+
+#define BOOST_TEST_MODULE alpha_shapes
+#include <boost/test/included/unit_test.hpp>
+#include <boost/system/error_code.hpp>
+#include <boost/chrono/thread_clock.hpp>
+// to construct a Delaunay_triangulation from a OFF file
+#include "gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h"
+#include "gudhi/Alpha_shapes.h"
+
+// to construct a simplex_tree from Delaunay_triangulation
+#include "gudhi/graph_simplicial_complex.h"
+#include "gudhi/Simplex_tree.h"
+
+#include <CGAL/Delaunay_triangulation.h>
+#include <CGAL/Epick_d.h>
+#include <CGAL/point_generators_d.h>
+#include <CGAL/algorithm.h>
+#include <CGAL/assertions.h>
+
+#include <iostream>
+#include <iterator>
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string>
+
+// Use dynamic_dimension_tag for the user to be able to set dimension
+typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > K;
+typedef CGAL::Delaunay_triangulation<K> T;
+// The triangulation uses the default instanciation of the 
+// TriangulationDataStructure template parameter
+
+BOOST_AUTO_TEST_CASE( OFF_file ) {
+  // ----------------------------------------------------------------------------
+  //
+  // Init of an alpha-shape from a OFF file
+  //
+  // ----------------------------------------------------------------------------
+  std::string off_file_name("S4_100.off");
+  std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl;
+  
+  Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_file(off_file_name, 4);
+
+  const int DIMENSION = 4;
+  std::cout << "alpha_shapes_from_file.dimension()=" << alpha_shapes_from_file.dimension() << std::endl;
+  BOOST_CHECK(alpha_shapes_from_file.dimension() == DIMENSION);
+  
+  const double FILTRATION = 0.0;
+  std::cout << "alpha_shapes_from_file.filtration()=" << alpha_shapes_from_file.filtration() << std::endl;
+  BOOST_CHECK(alpha_shapes_from_file.filtration() == FILTRATION);
+  
+  const int NUMBER_OF_VERTICES = 100;
+  std::cout << "alpha_shapes_from_file.num_vertices()=" << alpha_shapes_from_file.num_vertices() << std::endl;
+  BOOST_CHECK(alpha_shapes_from_file.num_vertices() == NUMBER_OF_VERTICES);
+  
+  const int NUMBER_OF_SIMPLICES = 6779;
+  std::cout << "alpha_shapes_from_file.num_simplices()=" << alpha_shapes_from_file.num_simplices() << std::endl;
+  BOOST_CHECK(alpha_shapes_from_file.num_simplices() == NUMBER_OF_SIMPLICES);
+
+}
+
+BOOST_AUTO_TEST_CASE( Delaunay_triangulation ) {
+  // ----------------------------------------------------------------------------
+  //
+  // Init of an alpha-shape from a Delauny triangulation
+  //
+  // ----------------------------------------------------------------------------
+  T dt(8);
+  std::string off_file_name("S8_10.off");
+  std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl;
+  
+  Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(off_file_name, dt, true, true);
+  std::cout << "off_reader.is_valid()=" << off_reader.is_valid() << std::endl;
+  BOOST_CHECK(off_reader.is_valid());
+
+  const int NUMBER_OF_VERTICES = 10;
+  std::cout << "dt.number_of_vertices()=" << dt.number_of_vertices() << std::endl;
+  BOOST_CHECK(dt.number_of_vertices() == NUMBER_OF_VERTICES);
+  
+  const int NUMBER_OF_FULL_CELLS = 30;
+  std::cout << "dt.number_of_full_cells()=" << dt.number_of_full_cells() << std::endl;
+  BOOST_CHECK(dt.number_of_full_cells() == NUMBER_OF_FULL_CELLS);
+  
+  const int NUMBER_OF_FINITE_FULL_CELLS = 6;
+  std::cout << "dt.number_of_finite_full_cells()=" << dt.number_of_finite_full_cells() << std::endl;
+  BOOST_CHECK(dt.number_of_finite_full_cells() == NUMBER_OF_FINITE_FULL_CELLS);
+
+  Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_dt(dt);
+
+  const int DIMENSION = 8;
+  std::cout << "alpha_shapes_from_dt.dimension()=" << alpha_shapes_from_dt.dimension() << std::endl;
+  BOOST_CHECK(alpha_shapes_from_dt.dimension() == DIMENSION);
+  
+  const double FILTRATION = 0.0;
+  std::cout << "alpha_shapes_from_dt.filtration()=" << alpha_shapes_from_dt.filtration() << std::endl;
+  BOOST_CHECK(alpha_shapes_from_dt.filtration() == FILTRATION);
+  
+  std::cout << "alpha_shapes_from_dt.num_vertices()=" << alpha_shapes_from_dt.num_vertices() << std::endl;
+  BOOST_CHECK(alpha_shapes_from_dt.num_vertices() == NUMBER_OF_VERTICES);
+  
+  const int NUMBER_OF_SIMPLICES = 997;
+  std::cout << "alpha_shapes_from_dt.num_simplices()=" << alpha_shapes_from_dt.num_simplices() << std::endl;
+  BOOST_CHECK(alpha_shapes_from_dt.num_simplices() == NUMBER_OF_SIMPLICES);
+}
+
diff --git a/src/Alpha_shapes/test/CMakeLists.txt b/src/Alpha_shapes/test/CMakeLists.txt
new file mode 100644
index 00000000..a48c1a8f
--- /dev/null
+++ b/src/Alpha_shapes/test/CMakeLists.txt
@@ -0,0 +1,31 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIAlphaShapesTest)
+
+# need CGAL 4.6
+# cmake -DCGAL_DIR=~/workspace/CGAL-4.6-beta1 ../../..
+if(CGAL_FOUND)
+  if (NOT CGAL_VERSION VERSION_LESS 4.6.0)
+    message(STATUS "CGAL version: ${CGAL_VERSION}.")
+
+    include( ${CGAL_USE_FILE} )
+ 
+    find_package(Eigen3 3.1.0)
+    if (EIGEN3_FOUND)
+      message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.")
+      include( ${EIGEN3_USE_FILE} )
+      include_directories (BEFORE "../../include")
+
+      add_definitions(-DDEBUG_TRACES)
+      add_executable ( AlphaShapesUnitTest Alpha_shapes_unit_test.cpp )
+      target_link_libraries(AlphaShapesUnitTest  ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+      add_test(AlphaShapesUnitTest ${CMAKE_CURRENT_BINARY_DIR}/AlphaShapesUnitTest)
+
+    else()
+      message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.")
+    endif()
+  else()
+    message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.")
+  endif ()
+endif()
+
+cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Alpha_shapes/include/gudhi")
diff --git a/src/Alpha_shapes/test/README b/src/Alpha_shapes/test/README
new file mode 100644
index 00000000..244a2b84
--- /dev/null
+++ b/src/Alpha_shapes/test/README
@@ -0,0 +1,12 @@
+To compile:
+***********
+
+cmake .
+make
+
+To launch with details:
+***********************
+
+./AlphaShapesUnitTest --report_level=detailed --log_level=all 
+
+   ==> echo $? returns 0 in case of success (non-zero otherwise)
diff --git a/src/Alpha_shapes/test/S4_100.off b/src/Alpha_shapes/test/S4_100.off
new file mode 100644
index 00000000..0a5dc58c
--- /dev/null
+++ b/src/Alpha_shapes/test/S4_100.off
@@ -0,0 +1,102 @@
+OFF
+100 0 0
+0.562921 -0.735261 -0.256472 0.277007
+-0.803733 -0.0527915 -0.315125 0.501918
+-0.24946 -0.354982 -0.410773 -0.801887
+0.916381 -0.0512295 0.371049 0.141223
+0.182222 0.940836 -0.171362 -0.228599
+-0.787145 -0.129213 0.568102 -0.202402
+0.0187866 -0.22093 -0.832882 -0.507095
+-0.4702 -0.533814 0.353776 0.607286
+-0.159798 -0.504771 0.263586 0.806346
+0.295546 0.162541 0.452931 0.825279
+0.242043 -0.107437 -0.913612 -0.308521
+0.875759 -0.113035 -0.469189 -0.0114505
+0.547877 -0.762247 -0.256972 -0.229729
+-0.172302 0.521057 0.412013 -0.727363
+-0.724729 -0.0574074 -0.0290602 0.686023
+0.700434 -0.102636 0.687285 0.162779
+-0.681386 0.0946893 0.610047 0.393178
+-0.847553 -0.357132 0.383743 0.0827718
+0.72297 -0.161631 -0.608517 0.284424
+0.757394 0.141549 0.196065 -0.606528
+0.78094 0.00901997 0.434536 0.448586
+0.14166 -0.619339 0.614589 -0.467582
+0.473105 -0.537832 -0.0103746 0.697711
+-0.208004 0.536218 0.818027 0.00605288
+0.743694 -0.628926 0.188072 0.126488
+-0.462228 -0.278147 0.35514 -0.76345
+-0.17361 0.249211 0.758567 -0.57648
+0.416958 -0.254924 -0.576373 -0.654946
+-0.590751 -0.286089 -0.424896 0.623402
+-0.639538 -0.739693 -0.203745 0.0482932
+0.0731787 0.132121 0.864022 0.480266
+-0.149644 -0.164724 -0.249746 -0.94239
+-0.348592 0.0120379 -0.928656 -0.126239
+0.395328 -0.54513 0.149976 -0.723917
+-0.974164 0.14707 0.157191 0.068302
+-0.166425 0.119943 -0.627304 -0.751269
+0.031947 -0.358518 -0.708301 -0.607251
+0.93781 -0.155368 -0.30951 0.0240237
+0.276094 -0.753155 -0.597088 -0.00387459
+-0.642876 -0.200551 -0.263517 -0.690687
+0.178711 0.604987 -0.262989 0.729993
+-0.520347 0.497922 -0.676144 0.155374
+-0.703999 0.500219 -0.484381 0.139789
+-0.131013 0.835735 0.506779 0.166004
+-0.536116 -0.566557 0.229226 0.582279
+-0.334105 0.158252 0.926091 0.0754059
+-0.0362677 0.296076 0.897108 0.325915
+-0.57486 0.798575 0.15324 -0.0912754
+0.498602 0.0186805 0.72824 0.469801
+-0.960329 0.0473356 0.261005 -0.0860505
+0.899134 -0.381392 -0.214508 0.00921711
+0.570576 0.567224 0.393019 -0.445237
+-0.761763 -0.614589 -0.0546476 -0.197513
+0.188584 0.289531 0.174031 0.922129
+-0.458506 -0.583876 0.639297 -0.2004
+0.785343 -0.21571 0.0794082 -0.574804
+0.0819036 0.65961 -0.247426 0.704973
+0.573125 0.49706 0.373026 0.534145
+-0.513286 -0.626226 0.208535 -0.548536
+0.460558 0.468686 0.507832 -0.55707
+0.716158 -0.488201 0.388209 -0.313164
+0.881074 0.152441 0.380128 -0.236589
+0.885793 0.0386389 0.161009 -0.433537
+-0.365162 0.298384 0.292846 0.831784
+0.364934 0.632269 -0.197205 -0.654346
+-0.31469 -0.429991 0.665304 -0.522923
+-0.734198 0.462914 -0.135691 -0.477756
+-0.422885 0.674444 -0.364143 -0.483419
+0.829218 -0.154622 -0.381147 0.378439
+0.887881 0.310479 -0.109528 0.321363
+-0.354398 -0.693974 0.456019 -0.429941
+-0.492045 -0.160008 0.044387 0.854587
+0.0595532 0.158421 0.412577 -0.895062
+-0.211441 0.491794 -0.153521 0.83058
+-0.33558 -0.504711 0.353831 -0.71236
+-0.735211 -0.197714 0.525626 0.379593
+0.465818 -0.424245 0.769469 -0.104627
+-0.641071 -0.286339 -0.704442 -0.103923
+-0.00446569 0.0249849 -0.194417 -0.980591
+-0.610081 -0.252448 0.176698 -0.729966
+-0.0859217 -0.154471 0.715027 0.676382
+0.091315 0.0723382 -0.855023 -0.505337
+0.165362 0.200983 -0.428242 -0.865373
+-0.587465 0.303019 -0.152442 0.734729
+0.454946 -0.319828 0.437063 -0.706902
+-0.384368 0.277509 0.879225 -0.0470385
+0.523335 -0.330233 -0.208592 0.757335
+0.895086 0.0448492 0.268089 -0.353466
+-0.0272491 -0.567336 -0.72254 -0.39411
+-0.0745014 -0.121818 -0.882466 0.448179
+0.382304 -0.240135 0.851109 -0.267941
+-0.418057 -0.852847 -0.3128 0.00606452
+-0.554046 0.304237 0.272381 -0.725453
+0.155115 -0.0894732 -0.245017 -0.952838
+0.114459 -0.130722 0.953669 0.245614
+0.0913002 -0.462466 0.244433 0.847374
+-0.198849 0.0785111 0.131441 -0.967997
+-0.303154 -0.686484 0.639333 0.167604
+0.521455 0.256835 -0.0584503 -0.811606
+-0.109787 0.870544 0.161523 0.451676
diff --git a/src/Alpha_shapes/test/S8_10.off b/src/Alpha_shapes/test/S8_10.off
new file mode 100644
index 00000000..1d67e10f
--- /dev/null
+++ b/src/Alpha_shapes/test/S8_10.off
@@ -0,0 +1,12 @@
+OFF
+10 0 0
+0.440036 -0.574754 -0.200485 0.216537 -0.501251 -0.0329236 -0.196529 0.313023
+-0.129367 -0.184089 -0.213021 -0.415848 0.783529 -0.0438025 0.317256 0.120749
+0.132429 0.683748 -0.124536 -0.166133 -0.540695 -0.0887576 0.390234 -0.139031
+0.0137399 -0.161581 -0.609142 -0.370872 -0.320669 -0.364053 0.24127 0.41416
+-0.115313 -0.36425 0.190208 0.581871 0.204605 0.112527 0.313562 0.571337
+0.168272 -0.0746917 -0.635156 -0.214488 0.629498 -0.0812499 -0.337255 -0.00823068
+0.369896 -0.514626 -0.173493 -0.1551 -0.127105 0.384377 0.303936 -0.536566
+-0.49013 -0.0388242 -0.0196532 0.463953 0.515962 -0.0756047 0.506276 0.119908
+-0.434258 0.060347 0.388793 0.250579 -0.653127 -0.275207 0.295714 0.0637842
+0.596172 -0.133284 -0.501793 0.234541 0.428452 0.0800735 0.110912 -0.343109
diff --git a/src/Bottleneck/concept/Persistence_diagram.h b/src/Bottleneck/concept/Persistence_diagram.h
new file mode 100644
index 00000000..eaaf8bc5
--- /dev/null
+++ b/src/Bottleneck/concept/Persistence_diagram.h
@@ -0,0 +1,7 @@
+typedef typename std::pair<double,double> Diagram_point;
+
+struct Persistence_Diagram
+{
+    const_iterator<Diagram_point> cbegin() const;
+    const_iterator<Diagram_point> cend() const;
+};
diff --git a/src/Bottleneck/example/CMakeLists.txt b/src/Bottleneck/example/CMakeLists.txt
new file mode 100644
index 00000000..2ff009c4
--- /dev/null
+++ b/src/Bottleneck/example/CMakeLists.txt
@@ -0,0 +1,5 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIBottleneckExample)
+
+add_executable ( RandomDiagrams random_diagrams.cpp )
+add_test(RandomDiagrams ${CMAKE_CURRENT_BINARY_DIR}/RandomDiagrams)
diff --git a/src/Bottleneck/example/random_diagrams.cpp b/src/Bottleneck/example/random_diagrams.cpp
new file mode 100644
index 00000000..71f152a6
--- /dev/null
+++ b/src/Bottleneck/example/random_diagrams.cpp
@@ -0,0 +1,40 @@
+/*    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):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (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/>.
+ */
+
+#include "gudhi/Graph_matching.h"
+#include <iostream>
+
+using namespace Gudhi::bottleneck;
+
+int main() {
+  int n = 100;
+  std::vector< std::pair<double, double> > v1, v2;
+  for (int i = 0; i < n; i++) {
+    int a = rand() % n;
+    v1.emplace_back(a, a + rand() % (n - a));
+    int b = rand() % n;
+    v2.emplace_back(b, b + rand() % (n - b));
+  }
+  // v1 and v2 are persistence diagrams containing each 100 randoms points.
+  double b = bottleneck_distance(v1, v2, 0);
+  std::cout << b << std::endl;
+}
diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h
new file mode 100644
index 00000000..ea47e1d5
--- /dev/null
+++ b/src/Bottleneck/include/gudhi/Graph_matching.h
@@ -0,0 +1,197 @@
+/*    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):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_
+
+#include <deque>
+#include <list>
+#include <vector>
+
+#include "gudhi/Layered_neighbors_finder.h"
+
+namespace Gudhi {
+
+namespace bottleneck {
+
+template<typename Persistence_diagram1, typename Persistence_diagram2>
+double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.);
+
+class Graph_matching {
+ public:
+  Graph_matching(const Persistence_diagrams_graph& g);
+  Graph_matching& operator=(const Graph_matching& m);
+  bool perfect() const;
+  bool multi_augment();
+  void set_r(double r);
+
+ private:
+  const Persistence_diagrams_graph& g;
+  double r;
+  std::vector<int> v_to_u;
+  std::list<int> unmatched_in_u;
+
+  Layered_neighbors_finder* layering() const;
+  bool augment(Layered_neighbors_finder* layered_nf, int u_start_index, int max_depth);
+  void update(std::deque<int>& path);
+};
+
+Graph_matching::Graph_matching(const Persistence_diagrams_graph& g)
+    : g(g), r(0), v_to_u(g.size()), unmatched_in_u() {
+  for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index)
+    unmatched_in_u.emplace_back(u_point_index);
+}
+
+Graph_matching& Graph_matching::operator=(const Graph_matching& m) {
+  r = m.r;
+  v_to_u = m.v_to_u;
+  unmatched_in_u = m.unmatched_in_u;
+  return *this;
+}
+
+inline bool Graph_matching::perfect() const {
+  return unmatched_in_u.empty();
+}
+
+inline bool Graph_matching::multi_augment() {
+  if (perfect())
+    return false;
+  Layered_neighbors_finder* layered_nf = layering();
+  double rn = sqrt(g.size());
+  int nblmax = layered_nf->vlayers_number()*2 + 1;
+  // verification of a necessary criterion
+  if ((unmatched_in_u.size() > rn && nblmax > rn) || nblmax == 0)
+    return false;
+  bool successful = false;
+  std::list<int>* tries = new std::list<int>(unmatched_in_u);
+  for (auto it = tries->cbegin(); it != tries->cend(); it++)
+    successful = successful || augment(layered_nf, *it, nblmax);
+  delete tries;
+  delete layered_nf;
+  return successful;
+}
+
+inline void Graph_matching::set_r(double r) {
+  this->r = r;
+}
+
+Layered_neighbors_finder* Graph_matching::layering() const {
+  bool end = false;
+  int layer = 0;
+  std::list<int> u_vertices(unmatched_in_u);
+  std::list<int> v_vertices;
+  Neighbors_finder nf(g, r);
+  Layered_neighbors_finder* layered_nf = new Layered_neighbors_finder(g, r);
+  for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index)
+    nf.add(v_point_index);
+  while (!u_vertices.empty()) {
+    for (auto it = u_vertices.cbegin(); it != u_vertices.cend(); ++it) {
+      std::list<int>* u_succ = nf.pull_all_near(*it);
+      for (auto it = u_succ->cbegin(); it != u_succ->cend(); ++it) {
+        layered_nf->add(*it, layer);
+        v_vertices.emplace_back(*it);
+      }
+      delete u_succ;
+    }
+    u_vertices.clear();
+    for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) {
+      if (v_to_u.at(*it) == null_point_index())
+        end = true;
+      else
+        u_vertices.emplace_back(v_to_u.at(*it));
+    }
+    if (end)
+      return layered_nf;
+    v_vertices.clear();
+    layer++;
+  }
+  return layered_nf;
+}
+
+bool Graph_matching::augment(Layered_neighbors_finder *layered_nf, int u_start_index, int max_depth) {
+  std::deque<int> path;
+  path.emplace_back(u_start_index);
+  // start is a point from U
+  do {
+    if (static_cast<int>(path.size()) > max_depth) {
+      path.pop_back();
+      path.pop_back();
+    }
+    if (path.empty())
+      return false;
+    int w = path.back();
+    path.emplace_back(layered_nf->pull_near(w, path.size() / 2));
+    while (path.back() == null_point_index()) {
+      path.pop_back();
+      path.pop_back();
+      if (path.empty())
+        return false;
+      path.pop_back();
+      path.emplace_back(layered_nf->pull_near(path.back(), path.size() / 2));
+    }
+    path.emplace_back(v_to_u.at(path.back()));
+  } while (path.back() != null_point_index());
+  path.pop_back();
+  update(path);
+  return true;
+}
+
+void Graph_matching::update(std::deque<int>& path) {
+  unmatched_in_u.remove(path.front());
+  for (auto it = path.cbegin(); it != path.cend(); ++it) {
+    int tmp = *it;
+    ++it;
+    v_to_u[*it] = tmp;
+  }
+}
+
+template<typename Persistence_diagram1, typename Persistence_diagram2>
+double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e) {
+  Persistence_diagrams_graph g(diag1, diag2, e);
+  std::vector<double>* sd = g.sorted_distances();
+  int idmin = 0;
+  int idmax = sd->size() - 1;
+  double alpha = pow(sd->size(), 0.25);
+  Graph_matching m(g);
+  Graph_matching biggest_unperfect = m;
+  while (idmin != idmax) {
+    int pas = static_cast<int>((idmax - idmin) / alpha);
+    m.set_r(sd->at(idmin + pas));
+    while (m.multi_augment()) {}
+    if (m.perfect()) {
+      idmax = idmin + pas;
+      m = biggest_unperfect;
+    } else {
+      biggest_unperfect = m;
+      idmin = idmin + pas + 1;
+    }
+  }
+  double b = sd->at(idmin);
+  delete sd;
+  return b;
+}
+
+}  // namespace bottleneck
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_
diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
new file mode 100644
index 00000000..de36e00b
--- /dev/null
+++ b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
@@ -0,0 +1,74 @@
+/*    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):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_
+
+#include <vector>
+
+#include "Neighbors_finder.h"
+
+// Layered_neighbors_finder is a data structure used to find if a query point from U has neighbors in V in a given
+// vlayer of the vlayered persistence diagrams graph. V's points have to be added manually using their index.
+// A neighbor returned is automatically removed.
+
+namespace Gudhi {
+
+namespace bottleneck {
+
+class Layered_neighbors_finder {
+ public:
+  Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r);
+  void add(int v_point_index, int vlayer);
+  int pull_near(int u_point_index, int vlayer);
+  int vlayers_number() const;
+
+ private:
+  const Persistence_diagrams_graph& g;
+  const double r;
+  std::vector<Neighbors_finder> neighbors_finder;
+};
+
+Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r) :
+    g(g), r(r), neighbors_finder() { }
+
+inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) {
+  for (int l = neighbors_finder.size(); l <= vlayer; l++)
+    neighbors_finder.emplace_back(Neighbors_finder(g, r));
+  neighbors_finder.at(vlayer).add(v_point_index);
+}
+
+inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) {
+  if (static_cast<int> (neighbors_finder.size()) <= vlayer)
+    return null_point_index();
+  return neighbors_finder.at(vlayer).pull_near(u_point_index);
+}
+
+inline int Layered_neighbors_finder::vlayers_number() const {
+  return neighbors_finder.size();
+}
+
+}  // namespace bottleneck
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_
diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h
new file mode 100644
index 00000000..98256571
--- /dev/null
+++ b/src/Bottleneck/include/gudhi/Neighbors_finder.h
@@ -0,0 +1,96 @@
+/*    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):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_
+
+#include <unordered_set>
+#include <list>
+
+#include "gudhi/Planar_neighbors_finder.h"
+
+namespace Gudhi {
+
+namespace bottleneck {
+
+// Neighbors_finder is a data structure used to find if a query point from U has neighbors in V
+// in the persistence diagrams graph.
+// V's points have to be added manually using their index. A neighbor returned is automatically removed.
+
+class Neighbors_finder {
+ public:
+  Neighbors_finder(const Persistence_diagrams_graph& g, double r);
+  void add(int v_point_index);
+  int pull_near(int u_point_index);
+  std::list<int>* pull_all_near(int u_point_index);
+
+ private:
+  const Persistence_diagrams_graph& g;
+  const double r;
+  Planar_neighbors_finder planar_neighbors_f;
+  std::unordered_set<int> projections_f;
+};
+
+Neighbors_finder::Neighbors_finder(const Persistence_diagrams_graph& g, double r) :
+    g(g), r(r), planar_neighbors_f(g, r), projections_f() { }
+
+inline void Neighbors_finder::add(int v_point_index) {
+  if (g.on_the_v_diagonal(v_point_index))
+    projections_f.emplace(v_point_index);
+  else
+    planar_neighbors_f.add(v_point_index);
+}
+
+inline int Neighbors_finder::pull_near(int u_point_index) {
+  int v_challenger = g.corresponding_point_in_v(u_point_index);
+  if (planar_neighbors_f.contains(v_challenger) && g.distance(u_point_index, v_challenger) < r) {
+    planar_neighbors_f.remove(v_challenger);
+    return v_challenger;
+  }
+  if (g.on_the_u_diagonal(u_point_index)) {
+    auto it = projections_f.cbegin();
+    if (it != projections_f.cend()) {
+      int tmp = *it;
+      projections_f.erase(it);
+      return tmp;
+    }
+  } else {
+    return planar_neighbors_f.pull_near(u_point_index);
+  }
+  return null_point_index();
+}
+
+inline std::list<int>* Neighbors_finder::pull_all_near(int u_point_index) {
+  std::list<int>* all_pull = planar_neighbors_f.pull_all_near(u_point_index);
+  int last_pull = pull_near(u_point_index);
+  while (last_pull != null_point_index()) {
+    all_pull->emplace_back(last_pull);
+    last_pull = pull_near(u_point_index);
+  }
+  return all_pull;
+}
+
+}  // namespace bottleneck
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_
diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
new file mode 100644
index 00000000..7e278209
--- /dev/null
+++ b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
@@ -0,0 +1,147 @@
+/*    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):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_
+
+#include <vector>
+#include <set>
+#include <cmath>
+#include <utility>  // for pair<>
+#include <algorithm>  // for max
+
+namespace Gudhi {
+
+namespace bottleneck {
+
+// Diagram_point is the type of the persistence diagram's points
+typedef typename std::pair<double, double> Diagram_point;
+
+// Return the used index for encoding none of the points
+int null_point_index();
+
+// Persistence_diagrams_graph is the interface beetwen any external representation of the two persistence diagrams and
+// the bottleneck distance computation. An interface is necessary to ensure basic functions complexity.
+
+class Persistence_diagrams_graph {
+ public:
+  // Persistence_diagram1 and 2 are the types of any externals representations of persistence diagrams.
+  // They have to have an iterator over points, which have to have fields first (for birth) and second (for death).
+  template<typename Persistence_diagram1, typename Persistence_diagram2>
+  Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.);
+  Persistence_diagrams_graph();
+  bool on_the_u_diagonal(int u_point_index) const;
+  bool on_the_v_diagonal(int v_point_index) const;
+  int corresponding_point_in_u(int v_point_index) const;
+  int corresponding_point_in_v(int u_point_index) const;
+  double distance(int u_point_index, int v_point_index) const;
+  int size() const;
+  std::vector<double>* sorted_distances();
+
+ private:
+  std::vector<Diagram_point> u;
+  std::vector<Diagram_point> v;
+  Diagram_point get_u_point(int u_point_index) const;
+  Diagram_point get_v_point(int v_point_index) const;
+};
+
+inline int null_point_index() {
+  return -1;
+}
+
+template<typename Persistence_diagram1, typename Persistence_diagram2>
+Persistence_diagrams_graph::Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e)
+    : u(), v() {
+  for (auto it = diag1.cbegin(); it != diag1.cend(); ++it)
+    if (it->second - it->first > e)
+      u.emplace_back(*it);
+  for (auto it = diag2.cbegin(); it != diag2.cend(); ++it)
+    if (it->second - it->first > e)
+      v.emplace_back(*it);
+  if (u.size() < v.size())
+    swap(u, v);
+}
+
+Persistence_diagrams_graph::Persistence_diagrams_graph::Persistence_diagrams_graph()
+    : u(), v() { }
+
+inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) const {
+  return u_point_index >= static_cast<int> (u.size());
+}
+
+inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) const {
+  return v_point_index >= static_cast<int> (v.size());
+}
+
+inline int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) const {
+  return on_the_v_diagonal(v_point_index) ?
+      v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size());
+}
+
+inline int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) const {
+  return on_the_u_diagonal(u_point_index) ?
+      u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size());
+}
+
+inline double Persistence_diagrams_graph::distance(int u_point_index, int v_point_index) const {
+  // could be optimized for the case where one point is the projection of the other
+  if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index))
+    return 0;
+  Diagram_point p_u = get_u_point(u_point_index);
+  Diagram_point p_v = get_v_point(v_point_index);
+  return std::max(std::fabs(p_u.first - p_v.first), std::fabs(p_u.second - p_v.second));
+}
+
+inline int Persistence_diagrams_graph::size() const {
+  return static_cast<int> (u.size() + v.size());
+}
+
+inline std::vector<double>* Persistence_diagrams_graph::sorted_distances() {
+  // could be optimized
+  std::set<double> sorted_distances;
+  for (int u_point_index = 0; u_point_index < size(); ++u_point_index)
+    for (int v_point_index = 0; v_point_index < size(); ++v_point_index)
+      sorted_distances.emplace(distance(u_point_index, v_point_index));
+  return new std::vector<double>(sorted_distances.cbegin(), sorted_distances.cend());
+}
+
+inline Diagram_point Persistence_diagrams_graph::get_u_point(int u_point_index) const {
+  if (!on_the_u_diagonal(u_point_index))
+    return u.at(u_point_index);
+  Diagram_point projector = v.at(corresponding_point_in_v(u_point_index));
+  double x = (projector.first + projector.second) / 2;
+  return Diagram_point(x, x);
+}
+
+inline Diagram_point Persistence_diagrams_graph::get_v_point(int v_point_index) const {
+  if (!on_the_v_diagonal(v_point_index))
+    return v.at(v_point_index);
+  Diagram_point projector = u.at(corresponding_point_in_u(v_point_index));
+  double x = (projector.first + projector.second) / 2;
+  return Diagram_point(x, x);
+}
+
+}  // namespace bottleneck
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_
diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
new file mode 100644
index 00000000..4af672e4
--- /dev/null
+++ b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
@@ -0,0 +1,119 @@
+/*    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):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_
+
+#include <list>
+#include <iostream>
+#include <set>
+
+#include "Persistence_diagrams_graph.h"
+
+namespace Gudhi {
+
+namespace bottleneck {
+
+// Planar_neighbors_finder is a data structure used to find if a query point from U has planar neighbors in V with the
+// planar distance.
+// V's points have to be added manually using their index. A neighbor returned is automatically removed but we can also
+// remove points manually using their index.
+
+class Abstract_planar_neighbors_finder {
+ public:
+  Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r);
+  virtual ~Abstract_planar_neighbors_finder() = 0;
+  virtual void add(int v_point_index) = 0;
+  virtual void remove(int v_point_index) = 0;
+  virtual bool contains(int v_point_index) const = 0;
+  virtual int pull_near(int u_point_index) = 0;
+  virtual std::list<int>* pull_all_near(int u_point_index);
+
+ protected:
+  const Persistence_diagrams_graph& g;
+  const double r;
+};
+
+
+// Naive_pnf is a nave implementation of Abstract_planar_neighbors_finder
+
+class Naive_pnf : public Abstract_planar_neighbors_finder {
+ public:
+  Naive_pnf(const Persistence_diagrams_graph& g, double r);
+  void add(int v_point_index);
+  void remove(int v_point_index);
+  bool contains(int v_point_index) const;
+  int pull_near(int u_point_index);
+
+ private:
+  std::set<int> candidates;
+};
+
+
+// Planar_neighbors_finder is the used Abstract_planar_neighbors_finder's implementation
+typedef Naive_pnf Planar_neighbors_finder;
+
+Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r) :
+    g(g), r(r) { }
+
+inline Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() { }
+
+inline std::list<int>* Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) {
+  std::list<int>* all_pull = new std::list<int>();
+  int last_pull = pull_near(u_point_index);
+  while (last_pull != null_point_index()) {
+    all_pull->emplace_back(last_pull);
+    last_pull = pull_near(u_point_index);
+  }
+  return all_pull;
+}
+
+Naive_pnf::Naive_pnf(const Persistence_diagrams_graph& g, double r) :
+    Abstract_planar_neighbors_finder(g, r), candidates() { }
+
+inline void Naive_pnf::add(int v_point_index) {
+  candidates.emplace(v_point_index);
+}
+
+inline void Naive_pnf::remove(int v_point_index) {
+  candidates.erase(v_point_index);
+}
+
+inline bool Naive_pnf::contains(int v_point_index) const {
+  return (candidates.count(v_point_index) > 0);
+}
+
+inline int Naive_pnf::pull_near(int u_point_index) {
+  for (auto it = candidates.begin(); it != candidates.end(); ++it)
+    if (g.distance(u_point_index, *it) <= r) {
+      int tmp = *it;
+      candidates.erase(it);
+      return tmp;
+    }
+  return null_point_index();
+}
+
+}  // namespace bottleneck
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_
diff --git a/src/Bottleneck/test/CMakeLists.txt b/src/Bottleneck/test/CMakeLists.txt
new file mode 100644
index 00000000..7044372e
--- /dev/null
+++ b/src/Bottleneck/test/CMakeLists.txt
@@ -0,0 +1,21 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIBottleneckUnitTest)
+
+if(NOT MSVC)
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --coverage")
+  set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} --coverage")	
+  set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} --coverage")
+endif()
+
+add_executable ( BottleneckUnitTest bottleneck_unit_test.cpp )
+target_link_libraries(BottleneckUnitTest ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+
+# Unitary tests
+add_test(BottleneckUnitTest ${CMAKE_CURRENT_BINARY_DIR}/BottleneckUnitTest)
+
+if (LCOV_PATH)
+  # Lcov code coverage of unitary test
+  add_test(src/Bottleneck/lcov/coverage.log ${CMAKE_SOURCE_DIR}/scripts/check_code_coverage.sh ${CMAKE_SOURCE_DIR}/src/Bottleneck)
+endif()
+
+cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Bottleneck/include/gudhi")
diff --git a/src/Bottleneck/test/README b/src/Bottleneck/test/README
new file mode 100644
index 00000000..0e7b8673
--- /dev/null
+++ b/src/Bottleneck/test/README
@@ -0,0 +1,12 @@
+To compile:
+***********
+
+cmake .
+make
+
+To launch with details:
+***********************
+
+./BottleneckUnitTest --report_level=detailed --log_level=all 
+
+   ==> echo $? returns 0 in case of success (non-zero otherwise)
diff --git a/src/Bottleneck/test/bottleneck_unit_test.cpp b/src/Bottleneck/test/bottleneck_unit_test.cpp
new file mode 100644
index 00000000..068b8690
--- /dev/null
+++ b/src/Bottleneck/test/bottleneck_unit_test.cpp
@@ -0,0 +1,26 @@
+#define BOOST_TEST_MODULE bottleneck test
+
+#include <boost/test/included/unit_test.hpp>
+
+#include "gudhi/Graph_matching.h"
+#include <iostream>
+
+using namespace Gudhi::bottleneck;
+
+BOOST_AUTO_TEST_CASE(random_diagrams) {
+  int n = 100;
+  // Random construction
+  std::vector< std::pair<double, double> > v1, v2;
+  for (int i = 0; i < n; i++) {
+    int a = rand() % n;
+    v1.emplace_back(a, a + rand() % (n - a));
+    int b = rand() % n;
+    v2.emplace_back(b, b + rand() % (n - b));
+  }
+  // v1 and v2 are persistence diagrams containing each 100 randoms points.
+  double b = bottleneck_distance(v1, v2, 0);
+  //
+  std::cout << b << std::endl;
+  const double EXPECTED_DISTANCE = 98.5;
+  BOOST_CHECK(b == EXPECTED_DISTANCE);
+}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 04c4369c..362ff7a8 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -38,5 +38,7 @@ else()
   add_subdirectory(example/Contraction)                   
   add_subdirectory(example/Hasse_complex)
   add_subdirectort(example/Witness_complex)
+  add_subdirectory(example/Alpha_shapes)
+  add_subdirectory(example/Bottleneck)
 
 endif()  
diff --git a/src/Hasse_complex/example/hasse_complex_from_simplex_tree.cpp b/src/Hasse_complex/example/hasse_complex_from_simplex_tree.cpp
index c7cff843..1de43ab7 100644
--- a/src/Hasse_complex/example/hasse_complex_from_simplex_tree.cpp
+++ b/src/Hasse_complex/example/hasse_complex_from_simplex_tree.cpp
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Vincent Rouvreau
  *
- *    Copyright (C) 2014  INRIA Sophia Antipolis-Méditerranée (France)
+ *    Copyright (C) 2014  INRIA Saclay (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
diff --git a/src/Simplex_tree/example/simple_simplex_tree.cpp b/src/Simplex_tree/example/simple_simplex_tree.cpp
index 4ecad752..bde224f1 100644
--- a/src/Simplex_tree/example/simple_simplex_tree.cpp
+++ b/src/Simplex_tree/example/simple_simplex_tree.cpp
@@ -31,273 +31,285 @@ typedef std::vector< Vertex_handle > typeVectorVertex;
 typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex;
 typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool;
 
-int main (int argc, char * const argv[])
-{
-	const Filtration_value FIRST_FILTRATION_VALUE = 0.1;
-	const Filtration_value SECOND_FILTRATION_VALUE = 0.2;
-	const Filtration_value THIRD_FILTRATION_VALUE = 0.3;
-	const Filtration_value FOURTH_FILTRATION_VALUE = 0.4;
-	Vertex_handle FIRST_VERTEX_HANDLE =  (Vertex_handle)0;
-	Vertex_handle SECOND_VERTEX_HANDLE = (Vertex_handle)1;
-	Vertex_handle THIRD_VERTEX_HANDLE = (Vertex_handle)2;
-	Vertex_handle FOURTH_VERTEX_HANDLE = (Vertex_handle)3;
-
-	// TEST OF INSERTION
-	std::cout << "********************************************************************" << std::endl;
-	std::cout << "EXAMPLE OF SIMPLE INSERTION" << std::endl;
-	//Construct the Simplex Tree
-	Simplex_tree<> simplexTree;
-
-	/* Simplex to be inserted:  */
-	/*    1                     */
-	/*    o                     */
-	/*   /X\                    */
-	/*  o---o---o               */
-	/*  2   0   3               */
-
-	// ++ FIRST
-	std::cout << "   * INSERT 0" << std::endl;
-	typeVectorVertex firstSimplexVector;
-	firstSimplexVector.push_back(FIRST_VERTEX_HANDLE);
-	typeSimplex firstSimplex = std::make_pair(firstSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
-	typePairSimplexBool returnValue =
-			simplexTree.insert_simplex ( firstSimplex.first, firstSimplex.second );
-
-	if (returnValue.second == true)
-	{
-		std::cout << "   + 0 INSERTED" << std::endl;
-		int nb_simplices = simplexTree.num_simplices() + 1;
-		simplexTree.set_num_simplices(nb_simplices);
-	}
-	else
-	{
-		std::cout << "   - 0 NOT INSERTED" << std::endl;
-	}
-
-	// ++ SECOND
-	std::cout << "   * INSERT 1" << std::endl;
-	typeVectorVertex secondSimplexVector;
-	secondSimplexVector.push_back(SECOND_VERTEX_HANDLE);
-	typeSimplex secondSimplex = std::make_pair(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
-	returnValue =
-			simplexTree.insert_simplex ( secondSimplex.first, secondSimplex.second );
-
-	if (returnValue.second == true)
-	{
-		std::cout << "   + 1 INSERTED" << std::endl;
-		int nb_simplices = simplexTree.num_simplices() + 1;
-		simplexTree.set_num_simplices(nb_simplices);
-	}
-	else
-	{
-		std::cout << "   - 1 NOT INSERTED" << std::endl;
-	}
-
-	// ++ THIRD
-	std::cout << "   * INSERT (0,1)" << std::endl;
-	typeVectorVertex thirdSimplexVector;
-	thirdSimplexVector.push_back(FIRST_VERTEX_HANDLE);
-	thirdSimplexVector.push_back(SECOND_VERTEX_HANDLE);
-	typeSimplex thirdSimplex = std::make_pair(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
-	returnValue =
-			simplexTree.insert_simplex ( thirdSimplex.first, thirdSimplex.second );
-
-	if (returnValue.second == true)
-	{
-		std::cout << "   + (0,1) INSERTED" << std::endl;
-		int nb_simplices = simplexTree.num_simplices() + 1;
-		simplexTree.set_num_simplices(nb_simplices);
-	}
-	else
-	{
-		std::cout << "   - (0,1) NOT INSERTED" << std::endl;
-	}
-
-	// ++ FOURTH
-	std::cout << "   * INSERT 2" << std::endl;
-	typeVectorVertex fourthSimplexVector;
-	fourthSimplexVector.push_back(THIRD_VERTEX_HANDLE);
-	typeSimplex fourthSimplex = std::make_pair(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
-	returnValue =
-			simplexTree.insert_simplex ( fourthSimplex.first, fourthSimplex.second );
-
-	if (returnValue.second == true)
-	{
-		std::cout << "   + 2 INSERTED" << std::endl;
-		int nb_simplices = simplexTree.num_simplices() + 1;
-		simplexTree.set_num_simplices(nb_simplices);
-	}
-	else
-	{
-		std::cout << "   - 2 NOT INSERTED" << std::endl;
-	}
-
-	// ++ FIFTH
-	std::cout << "   * INSERT (2,0)" << std::endl;
-	typeVectorVertex fifthSimplexVector;
-	fifthSimplexVector.push_back(THIRD_VERTEX_HANDLE);
-	fifthSimplexVector.push_back(FIRST_VERTEX_HANDLE);
-	typeSimplex fifthSimplex = std::make_pair(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
-	returnValue =
-			simplexTree.insert_simplex ( fifthSimplex.first, fifthSimplex.second );
-
-	if (returnValue.second == true)
-	{
-		std::cout << "   + (2,0) INSERTED" << std::endl;
-		int nb_simplices = simplexTree.num_simplices() + 1;
-		simplexTree.set_num_simplices(nb_simplices);
-	}
-	else
-	{
-		std::cout << "   - (2,0) NOT INSERTED" << std::endl;
-	}
-
-	// ++ SIXTH
-	std::cout << "   * INSERT (2,1)" << std::endl;
-	typeVectorVertex sixthSimplexVector;
-	sixthSimplexVector.push_back(THIRD_VERTEX_HANDLE);
-	sixthSimplexVector.push_back(SECOND_VERTEX_HANDLE);
-	typeSimplex sixthSimplex = std::make_pair(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
-	returnValue =
-			simplexTree.insert_simplex ( sixthSimplex.first, sixthSimplex.second );
-
-	if (returnValue.second == true)
-	{
-		std::cout << "   + (2,1) INSERTED" << std::endl;
-		int nb_simplices = simplexTree.num_simplices() + 1;
-		simplexTree.set_num_simplices(nb_simplices);
-	}
-	else
-	{
-		std::cout << "   - (2,1) NOT INSERTED" << std::endl;
-	}
-
-	// ++ SEVENTH
-	std::cout << "   * INSERT (2,1,0)" << std::endl;
-	typeVectorVertex seventhSimplexVector;
-	seventhSimplexVector.push_back(THIRD_VERTEX_HANDLE);
-	seventhSimplexVector.push_back(SECOND_VERTEX_HANDLE);
-	seventhSimplexVector.push_back(FIRST_VERTEX_HANDLE);
-	typeSimplex seventhSimplex = std::make_pair(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE));
-	returnValue =
-			simplexTree.insert_simplex ( seventhSimplex.first, seventhSimplex.second );
-
-	if (returnValue.second == true)
-	{
-		std::cout << "   + (2,1,0) INSERTED" << std::endl;
-		int nb_simplices = simplexTree.num_simplices() + 1;
-		simplexTree.set_num_simplices(nb_simplices);
-	}
-	else
-	{
-		std::cout << "   - (2,1,0) NOT INSERTED" << std::endl;
-	}
-
-	// ++ EIGHTH
-	std::cout << "   * INSERT 3" << std::endl;
-	typeVectorVertex eighthSimplexVector;
-	eighthSimplexVector.push_back(FOURTH_VERTEX_HANDLE);
-	typeSimplex eighthSimplex = std::make_pair(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
-	returnValue =
-			simplexTree.insert_simplex ( eighthSimplex.first, eighthSimplex.second );
-
-	if (returnValue.second == true)
-	{
-		std::cout << "   + 3 INSERTED" << std::endl;
-		int nb_simplices = simplexTree.num_simplices() + 1;
-		simplexTree.set_num_simplices(nb_simplices);
-	}
-	else
-	{
-		std::cout << "   - 3 NOT INSERTED" << std::endl;
-	}
-
-	// ++ NINETH
-	std::cout << "   * INSERT (3,0)" << std::endl;
-	typeVectorVertex ninethSimplexVector;
-	ninethSimplexVector.push_back(FOURTH_VERTEX_HANDLE);
-	ninethSimplexVector.push_back(FIRST_VERTEX_HANDLE);
-	typeSimplex ninethSimplex = std::make_pair(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
-	returnValue =
-			simplexTree.insert_simplex ( ninethSimplex.first, ninethSimplex.second );
-
-	if (returnValue.second == true)
-	{
-		std::cout << "   + (3,0) INSERTED" << std::endl;
-		int nb_simplices = simplexTree.num_simplices() + 1;
-		simplexTree.set_num_simplices(nb_simplices);
-	}
-	else
-	{
-		std::cout << "   - (3,0) NOT INSERTED" << std::endl;
-	}
-
-	// ++ TENTH
-	std::cout << "   * INSERT 0 (already inserted)" << std::endl;
-	typeVectorVertex tenthSimplexVector;
-	tenthSimplexVector.push_back(FIRST_VERTEX_HANDLE);
-	typeSimplex tenthSimplex = std::make_pair(tenthSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); // With a different filtration value
-	returnValue =
-			simplexTree.insert_simplex ( tenthSimplex.first, tenthSimplex.second );
-
-	if (returnValue.second == true)
-	{
-		std::cout << "   + 0 INSERTED" << std::endl;
-		int nb_simplices = simplexTree.num_simplices() + 1;
-		simplexTree.set_num_simplices(nb_simplices);
-	}
-	else
-	{
-		std::cout << "   - 0 NOT INSERTED" << std::endl;
-	}
-
-	// ++ ELEVENTH
-	std::cout << "   * INSERT (2,1,0) (already inserted)" << std::endl;
-	typeVectorVertex eleventhSimplexVector;
-	eleventhSimplexVector.push_back(THIRD_VERTEX_HANDLE);
-	eleventhSimplexVector.push_back(SECOND_VERTEX_HANDLE);
-	eleventhSimplexVector.push_back(FIRST_VERTEX_HANDLE);
-	typeSimplex eleventhSimplex = std::make_pair(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE));
-	returnValue =
-			simplexTree.insert_simplex ( eleventhSimplex.first, eleventhSimplex.second );
-
-	if (returnValue.second == true)
-	{
-		std::cout << "   + (2,1,0) INSERTED" << std::endl;
-		int nb_simplices = simplexTree.num_simplices() + 1;
-		simplexTree.set_num_simplices(nb_simplices);
-	}
-	else
-	{
-		std::cout << "   - (2,1,0) NOT INSERTED" << std::endl;
-	}
-
-	// ++ GENERAL VARIABLE SET
-	simplexTree.set_filtration(FOURTH_FILTRATION_VALUE); // Max filtration value
-	simplexTree.set_dimension(2); // Max dimension = 2 -> (2,1,0)
-
-	std::cout << "********************************************************************" << std::endl;
-	// Display the Simplex_tree - Can not be done in the middle of 2 inserts
-	std::cout << "* The complex contains " << simplexTree.num_simplices() << " simplices" << std::endl;
-	std::cout << "   - dimension " << simplexTree.dimension() << "   - filtration " << simplexTree.filtration() << std::endl;
-	std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:" << std::endl;
-	for( auto f_simplex : simplexTree.filtration_simplex_range() )
-	{
-		std::cout << "   " << "[" << simplexTree.filtration(f_simplex) << "] ";
-		for( auto vertex : simplexTree.simplex_vertex_range(f_simplex) )
-		{
-			std::cout << (int)vertex << " ";
-		}
-		std::cout << std::endl;
-	}
-	//   [0.1] 0
-	//   [0.1] 1
-	//   [0.1] 2
-	//   [0.1] 3
-	//   [0.2] 1 0
-	//   [0.2] 2 0
-	//   [0.2] 2 1
-	//   [0.2] 3 0
-	//   [0.3] 2 1 0
-	return 0;
+int main(int argc, char * const argv[]) {
+  const Filtration_value FIRST_FILTRATION_VALUE = 0.1;
+  const Filtration_value SECOND_FILTRATION_VALUE = 0.2;
+  const Filtration_value THIRD_FILTRATION_VALUE = 0.3;
+  const Filtration_value FOURTH_FILTRATION_VALUE = 0.4;
+  Vertex_handle FIRST_VERTEX_HANDLE = (Vertex_handle) 0;
+  Vertex_handle SECOND_VERTEX_HANDLE = (Vertex_handle) 1;
+  Vertex_handle THIRD_VERTEX_HANDLE = (Vertex_handle) 2;
+  Vertex_handle FOURTH_VERTEX_HANDLE = (Vertex_handle) 3;
+
+  // TEST OF INSERTION
+  std::cout << "********************************************************************" << std::endl;
+  std::cout << "EXAMPLE OF SIMPLE INSERTION" << std::endl;
+  //Construct the Simplex Tree
+  Simplex_tree<> simplexTree;
+
+  /* Simplex to be inserted:  */
+  /*    1                     */
+  /*    o                     */
+  /*   /X\                    */
+  /*  o---o---o               */
+  /*  2   0   3               */
+
+  // ++ FIRST
+  std::cout << "   * INSERT 0" << std::endl;
+  typeVectorVertex firstSimplexVector;
+  firstSimplexVector.push_back(FIRST_VERTEX_HANDLE);
+  typeSimplex firstSimplex = std::make_pair(firstSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
+  typePairSimplexBool returnValue =
+      simplexTree.insert_simplex(firstSimplex.first, firstSimplex.second);
+
+  if (returnValue.second == true) {
+    std::cout << "   + 0 INSERTED" << std::endl;
+    int nb_simplices = simplexTree.num_simplices() + 1;
+    simplexTree.set_num_simplices(nb_simplices);
+  } else {
+    std::cout << "   - 0 NOT INSERTED" << std::endl;
+  }
+
+  // ++ SECOND
+  std::cout << "   * INSERT 1" << std::endl;
+  typeVectorVertex secondSimplexVector;
+  secondSimplexVector.push_back(SECOND_VERTEX_HANDLE);
+  typeSimplex secondSimplex = std::make_pair(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
+  returnValue =
+      simplexTree.insert_simplex(secondSimplex.first, secondSimplex.second);
+
+  if (returnValue.second == true) {
+    std::cout << "   + 1 INSERTED" << std::endl;
+    int nb_simplices = simplexTree.num_simplices() + 1;
+    simplexTree.set_num_simplices(nb_simplices);
+  } else {
+    std::cout << "   - 1 NOT INSERTED" << std::endl;
+  }
+
+  // ++ THIRD
+  std::cout << "   * INSERT (0,1)" << std::endl;
+  typeVectorVertex thirdSimplexVector;
+  thirdSimplexVector.push_back(FIRST_VERTEX_HANDLE);
+  thirdSimplexVector.push_back(SECOND_VERTEX_HANDLE);
+  typeSimplex thirdSimplex = std::make_pair(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  returnValue =
+      simplexTree.insert_simplex(thirdSimplex.first, thirdSimplex.second);
+
+  if (returnValue.second == true) {
+    std::cout << "   + (0,1) INSERTED" << std::endl;
+    int nb_simplices = simplexTree.num_simplices() + 1;
+    simplexTree.set_num_simplices(nb_simplices);
+  } else {
+    std::cout << "   - (0,1) NOT INSERTED" << std::endl;
+  }
+
+  // ++ FOURTH
+  std::cout << "   * INSERT 2" << std::endl;
+  typeVectorVertex fourthSimplexVector;
+  fourthSimplexVector.push_back(THIRD_VERTEX_HANDLE);
+  typeSimplex fourthSimplex = std::make_pair(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
+  returnValue =
+      simplexTree.insert_simplex(fourthSimplex.first, fourthSimplex.second);
+
+  if (returnValue.second == true) {
+    std::cout << "   + 2 INSERTED" << std::endl;
+    int nb_simplices = simplexTree.num_simplices() + 1;
+    simplexTree.set_num_simplices(nb_simplices);
+  } else {
+    std::cout << "   - 2 NOT INSERTED" << std::endl;
+  }
+
+  // ++ FIFTH
+  std::cout << "   * INSERT (2,0)" << std::endl;
+  typeVectorVertex fifthSimplexVector;
+  fifthSimplexVector.push_back(THIRD_VERTEX_HANDLE);
+  fifthSimplexVector.push_back(FIRST_VERTEX_HANDLE);
+  typeSimplex fifthSimplex = std::make_pair(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  returnValue =
+      simplexTree.insert_simplex(fifthSimplex.first, fifthSimplex.second);
+
+  if (returnValue.second == true) {
+    std::cout << "   + (2,0) INSERTED" << std::endl;
+    int nb_simplices = simplexTree.num_simplices() + 1;
+    simplexTree.set_num_simplices(nb_simplices);
+  } else {
+    std::cout << "   - (2,0) NOT INSERTED" << std::endl;
+  }
+
+  // ++ SIXTH
+  std::cout << "   * INSERT (2,1)" << std::endl;
+  typeVectorVertex sixthSimplexVector;
+  sixthSimplexVector.push_back(THIRD_VERTEX_HANDLE);
+  sixthSimplexVector.push_back(SECOND_VERTEX_HANDLE);
+  typeSimplex sixthSimplex = std::make_pair(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  returnValue =
+      simplexTree.insert_simplex(sixthSimplex.first, sixthSimplex.second);
+
+  if (returnValue.second == true) {
+    std::cout << "   + (2,1) INSERTED" << std::endl;
+    int nb_simplices = simplexTree.num_simplices() + 1;
+    simplexTree.set_num_simplices(nb_simplices);
+  } else {
+    std::cout << "   - (2,1) NOT INSERTED" << std::endl;
+  }
+
+  // ++ SEVENTH
+  std::cout << "   * INSERT (2,1,0)" << std::endl;
+  typeVectorVertex seventhSimplexVector;
+  seventhSimplexVector.push_back(THIRD_VERTEX_HANDLE);
+  seventhSimplexVector.push_back(SECOND_VERTEX_HANDLE);
+  seventhSimplexVector.push_back(FIRST_VERTEX_HANDLE);
+  typeSimplex seventhSimplex = std::make_pair(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE));
+  returnValue =
+      simplexTree.insert_simplex(seventhSimplex.first, seventhSimplex.second);
+
+  if (returnValue.second == true) {
+    std::cout << "   + (2,1,0) INSERTED" << std::endl;
+    int nb_simplices = simplexTree.num_simplices() + 1;
+    simplexTree.set_num_simplices(nb_simplices);
+  } else {
+    std::cout << "   - (2,1,0) NOT INSERTED" << std::endl;
+  }
+
+  // ++ EIGHTH
+  std::cout << "   * INSERT 3" << std::endl;
+  typeVectorVertex eighthSimplexVector;
+  eighthSimplexVector.push_back(FOURTH_VERTEX_HANDLE);
+  typeSimplex eighthSimplex = std::make_pair(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
+  returnValue =
+      simplexTree.insert_simplex(eighthSimplex.first, eighthSimplex.second);
+
+  if (returnValue.second == true) {
+    std::cout << "   + 3 INSERTED" << std::endl;
+    int nb_simplices = simplexTree.num_simplices() + 1;
+    simplexTree.set_num_simplices(nb_simplices);
+  } else {
+    std::cout << "   - 3 NOT INSERTED" << std::endl;
+  }
+
+  // ++ NINETH
+  std::cout << "   * INSERT (3,0)" << std::endl;
+  typeVectorVertex ninethSimplexVector;
+  ninethSimplexVector.push_back(FOURTH_VERTEX_HANDLE);
+  ninethSimplexVector.push_back(FIRST_VERTEX_HANDLE);
+  typeSimplex ninethSimplex = std::make_pair(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  returnValue =
+      simplexTree.insert_simplex(ninethSimplex.first, ninethSimplex.second);
+
+  if (returnValue.second == true) {
+    std::cout << "   + (3,0) INSERTED" << std::endl;
+    int nb_simplices = simplexTree.num_simplices() + 1;
+    simplexTree.set_num_simplices(nb_simplices);
+  } else {
+    std::cout << "   - (3,0) NOT INSERTED" << std::endl;
+  }
+
+  // ++ TENTH
+  std::cout << "   * INSERT 0 (already inserted)" << std::endl;
+  typeVectorVertex tenthSimplexVector;
+  tenthSimplexVector.push_back(FIRST_VERTEX_HANDLE);
+  typeSimplex tenthSimplex = std::make_pair(tenthSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE)); // With a different filtration value
+  returnValue =
+      simplexTree.insert_simplex(tenthSimplex.first, tenthSimplex.second);
+
+  if (returnValue.second == true) {
+    std::cout << "   + 0 INSERTED" << std::endl;
+    int nb_simplices = simplexTree.num_simplices() + 1;
+    simplexTree.set_num_simplices(nb_simplices);
+  } else {
+    std::cout << "   - 0 NOT INSERTED" << std::endl;
+  }
+
+  // ++ ELEVENTH
+  std::cout << "   * INSERT (2,1,0) (already inserted)" << std::endl;
+  typeVectorVertex eleventhSimplexVector;
+  eleventhSimplexVector.push_back(THIRD_VERTEX_HANDLE);
+  eleventhSimplexVector.push_back(SECOND_VERTEX_HANDLE);
+  eleventhSimplexVector.push_back(FIRST_VERTEX_HANDLE);
+  typeSimplex eleventhSimplex = std::make_pair(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE));
+  returnValue =
+      simplexTree.insert_simplex(eleventhSimplex.first, eleventhSimplex.second);
+
+  if (returnValue.second == true) {
+    std::cout << "   + (2,1,0) INSERTED" << std::endl;
+    int nb_simplices = simplexTree.num_simplices() + 1;
+    simplexTree.set_num_simplices(nb_simplices);
+  } else {
+    std::cout << "   - (2,1,0) NOT INSERTED" << std::endl;
+  }
+
+  // ++ GENERAL VARIABLE SET
+  simplexTree.set_filtration(FOURTH_FILTRATION_VALUE); // Max filtration value
+  simplexTree.set_dimension(2); // Max dimension = 2 -> (2,1,0)
+
+  std::cout << "********************************************************************" << std::endl;
+  // Display the Simplex_tree - Can not be done in the middle of 2 inserts
+  std::cout << "* The complex contains " << simplexTree.num_simplices() << " simplices" << std::endl;
+  std::cout << "   - dimension " << simplexTree.dimension() << "   - filtration " << simplexTree.filtration() << std::endl;
+  std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:" << std::endl;
+  for (auto f_simplex : simplexTree.filtration_simplex_range()) {
+    std::cout << "   " << "[" << simplexTree.filtration(f_simplex) << "] ";
+    for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) {
+      std::cout << (int) vertex << " ";
+    }
+    std::cout << std::endl;
+  }
+  //   [0.1] 0
+  //   [0.1] 1
+  //   [0.1] 2
+  //   [0.1] 3
+  //   [0.2] 1 0
+  //   [0.2] 2 0
+  //   [0.2] 2 1
+  //   [0.2] 3 0
+  //   [0.3] 2 1 0
+
+  // ------------------------------------------------------------------------------------------------------------------
+  // Find in the simplex_tree
+  // ------------------------------------------------------------------------------------------------------------------
+  Simplex_tree<>::Simplex_handle simplexFound = simplexTree.find(secondSimplexVector);
+  std::cout << "**************IS THE SIMPLEX {1} IN THE SIMPLEX TREE ?\n";
+  if (simplexFound != simplexTree.null_simplex())
+    std::cout << "***+ YES IT IS!\n";
+  else
+    std::cout << "***- NO IT ISN'T\n";
+
+  Vertex_handle UNKNOWN_VERTEX_HANDLE = (Vertex_handle) 15;
+  typeVectorVertex unknownSimplexVector;
+  unknownSimplexVector.push_back(UNKNOWN_VERTEX_HANDLE);
+  simplexFound = simplexTree.find(unknownSimplexVector);
+  std::cout << "**************IS THE SIMPLEX {15} IN THE SIMPLEX TREE ?\n";
+  if (simplexFound != simplexTree.null_simplex())
+    std::cout << "***+ YES IT IS!\n";
+  else
+    std::cout << "***- NO IT ISN'T\n";
+
+  simplexFound = simplexTree.find(fifthSimplexVector);
+  std::cout << "**************IS THE SIMPLEX {2,0} IN THE SIMPLEX TREE ?\n";
+  if (simplexFound != simplexTree.null_simplex())
+    std::cout << "***+ YES IT IS!\n";
+  else
+    std::cout << "***- NO IT ISN'T\n";
+
+  typeVectorVertex otherSimplexVector;
+  otherSimplexVector.push_back(UNKNOWN_VERTEX_HANDLE);
+  otherSimplexVector.push_back(SECOND_VERTEX_HANDLE);
+  simplexFound = simplexTree.find(otherSimplexVector);
+  std::cout << "**************IS THE SIMPLEX {15,1} IN THE SIMPLEX TREE ?\n";
+  if (simplexFound != simplexTree.null_simplex())
+    std::cout << "***+ YES IT IS!\n";
+  else
+    std::cout << "***- NO IT ISN'T\n";
+
+  typeVectorVertex invSimplexVector;
+  invSimplexVector.push_back(SECOND_VERTEX_HANDLE);
+  invSimplexVector.push_back(THIRD_VERTEX_HANDLE);
+  invSimplexVector.push_back(FIRST_VERTEX_HANDLE);
+  simplexFound = simplexTree.find(invSimplexVector);
+  std::cout << "**************IS THE SIMPLEX {1,2,0} IN THE SIMPLEX TREE ?\n";
+  if (simplexFound != simplexTree.null_simplex())
+    std::cout << "***+ YES IT IS!\n";
+  else
+    std::cout << "***- NO IT ISN'T\n";
+  return 0;
 }
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index e52fe1ae..f95679cb 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -397,7 +397,7 @@ class Simplex_tree {
    * <CODE>Vertex_handle</CODE>.
    */
   template<class RandomAccessVertexRange>
-  Simplex_handle find(const RandomAccessVertexRange & s) {
+  Simplex_handle find(RandomAccessVertexRange & s) {
     if (s.begin() == s.end())
       std::cerr << "Empty simplex \n";
 
diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
index eaec4881..c0cfced1 100644
--- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp
+++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
@@ -512,6 +512,67 @@ BOOST_AUTO_TEST_CASE( NSimplexAndSubfaces_tree_insertion )
   test_simplex_tree_contains(st,simplexPair4,2);  // (1,0) is in position 2
   test_simplex_tree_contains(st,simplexPair5,14);  // (3,4,5) is in position 14
   test_simplex_tree_contains(st,simplexPair6,26);  // (7,6,1,0) is in position 26
+  
+  // ------------------------------------------------------------------------------------------------------------------
+  // Find in the simplex_tree
+  // ------------------------------------------------------------------------------------------------------------------
+  typeVectorVertex simpleSimplexVector;
+  simpleSimplexVector.push_back(SECOND_VERTEX_HANDLE);
+  Simplex_tree<>::Simplex_handle simplexFound = st.find(simpleSimplexVector);
+  std::cout << "**************IS THE SIMPLEX {1} IN THE SIMPLEX TREE ?\n";
+  if (simplexFound != st.null_simplex())
+    std::cout << "***+ YES IT IS!\n";
+  else
+    std::cout << "***- NO IT ISN'T\n";
+  // Check it is found
+  BOOST_CHECK(simplexFound != st.null_simplex());
+
+  Vertex_handle UNKNOWN_VERTEX_HANDLE = (Vertex_handle) 15;
+  typeVectorVertex unknownSimplexVector;
+  unknownSimplexVector.push_back(UNKNOWN_VERTEX_HANDLE);
+  simplexFound = st.find(unknownSimplexVector);
+  std::cout << "**************IS THE SIMPLEX {15} IN THE SIMPLEX TREE ?\n";
+  if (simplexFound != st.null_simplex())
+    std::cout << "***+ YES IT IS!\n";
+  else
+    std::cout << "***- NO IT ISN'T\n";
+  // Check it is NOT found
+  BOOST_CHECK(simplexFound == st.null_simplex());
+
+  simplexFound = st.find(SimplexVector6);
+  std::cout << "**************IS THE SIMPLEX {0,1,6,7} IN THE SIMPLEX TREE ?\n";
+  if (simplexFound != st.null_simplex())
+    std::cout << "***+ YES IT IS!\n";
+  else
+    std::cout << "***- NO IT ISN'T\n";
+  // Check it is found
+  BOOST_CHECK(simplexFound != st.null_simplex());
+      
+  typeVectorVertex otherSimplexVector;
+  otherSimplexVector.push_back(UNKNOWN_VERTEX_HANDLE);
+  otherSimplexVector.push_back(SECOND_VERTEX_HANDLE);
+  simplexFound = st.find(otherSimplexVector);
+  std::cout << "**************IS THE SIMPLEX {15,1} IN THE SIMPLEX TREE ?\n";
+  if (simplexFound != st.null_simplex())
+    std::cout << "***+ YES IT IS!\n";
+  else
+    std::cout << "***- NO IT ISN'T\n";
+  // Check it is NOT found
+  BOOST_CHECK(simplexFound == st.null_simplex());
+
+  typeVectorVertex invSimplexVector;
+  invSimplexVector.push_back(SECOND_VERTEX_HANDLE);
+  invSimplexVector.push_back(THIRD_VERTEX_HANDLE);
+  invSimplexVector.push_back(FIRST_VERTEX_HANDLE);
+  simplexFound = st.find(invSimplexVector);
+  std::cout << "**************IS THE SIMPLEX {1,2,0} IN THE SIMPLEX TREE ?\n";
+  if (simplexFound != st.null_simplex())
+    std::cout << "***+ YES IT IS!\n";
+  else
+    std::cout << "***- NO IT ISN'T\n";
+  // Check it is found
+  BOOST_CHECK(simplexFound != st.null_simplex());
+  
   // Display the Simplex_tree - Can not be done in the middle of 2 inserts
   std::cout << "The complex contains " << st.num_simplices() << " simplices" << std::endl;
   std::cout << "   - dimension " << st.dimension() << "   - filtration " << st.filtration() << std::endl;
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h
index c98b0b45..aaaab8b0 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h
@@ -24,6 +24,7 @@
 
 #include <string>
 #include <vector>
+#include <map>
 
 #include "gudhi/Off_reader.h"
 
@@ -36,38 +37,36 @@ namespace skbl {
  */
 template<typename Complex>
 class Skeleton_blocker_off_flag_visitor_reader {
-	Complex& complex_;
-	typedef typename Complex::Vertex_handle Vertex_handle;
-	typedef typename Complex::Point Point;
-
-	const bool load_only_points_;
-
-public:
-	explicit Skeleton_blocker_off_flag_visitor_reader(Complex& complex, bool load_only_points = false):
-	complex_(complex),
-	load_only_points_(load_only_points) {}
-
-
-	void init(int dim, int num_vertices, int num_faces, int num_edges) {
-		// todo do an assert to check that this number are correctly read
-		// todo reserve size for vector points
-	}
-
-
-	void point(const std::vector<double>& point) {
-		complex_.add_vertex(Point(point.begin(),point.end()));
-	}
-
-	void maximal_face(const std::vector<int>& face) {
-		if (!load_only_points_) {
-			for (size_t i = 0; i < face.size(); ++i)
-				for (size_t j = i+1; j < face.size(); ++j) {
-					complex_.add_edge(Vertex_handle(face[i]), Vertex_handle(face[j]));
-				}
-		}
-	}
-
-	void done() {}
+  Complex& complex_;
+  typedef typename Complex::Vertex_handle Vertex_handle;
+  typedef typename Complex::Point Point;
+
+  const bool load_only_points_;
+
+ public:
+  explicit Skeleton_blocker_off_flag_visitor_reader(Complex& complex, bool load_only_points = false) :
+      complex_(complex),
+      load_only_points_(load_only_points) { }
+
+  void init(int dim, int num_vertices, int num_faces, int num_edges) {
+    // todo do an assert to check that this number are correctly read
+    // todo reserve size for vector points
+  }
+
+  void point(const std::vector<double>& point) {
+    complex_.add_vertex(Point(point.begin(), point.end()));
+  }
+
+  void maximal_face(const std::vector<int>& face) {
+    if (!load_only_points_) {
+      for (size_t i = 0; i < face.size(); ++i)
+        for (size_t j = i + 1; j < face.size(); ++j) {
+          complex_.add_edge(Vertex_handle(face[i]), Vertex_handle(face[j]));
+        }
+    }
+  }
+
+  void done() { }
 };
 
 /**
@@ -75,137 +74,127 @@ public:
  */
 template<typename Complex>
 class Skeleton_blocker_off_visitor_reader {
-	Complex& complex_;
-	typedef typename Complex::Vertex_handle Vertex_handle;
-	typedef typename Complex::Simplex_handle Simplex_handle;
-	typedef typename Complex::Point Point;
-
-	const bool load_only_points_;
-	std::vector<Point> points_;
-	std::vector<Simplex_handle> maximal_faces_;
-
-public:
-	explicit Skeleton_blocker_off_visitor_reader(Complex& complex, bool load_only_points = false):
-	complex_(complex),
-	load_only_points_(load_only_points) {}
-
-
-	void init(int dim, int num_vertices, int num_faces, int num_edges){
-		maximal_faces_.reserve(num_faces);
-		points_.reserve(num_vertices);
-	}
-
-
-	void point(const std::vector<double>& point) {
-		points_.emplace_back(point.begin(),point.end());
-	}
-
-	void maximal_face(const std::vector<int>& face) {
-		if (!load_only_points_) {
-			Simplex_handle s;
-			for(auto x : face)
-				s.add_vertex(Vertex_handle(x));
-			maximal_faces_.emplace_back(s);
-		}
-	}
-
-	void done() {
-		complex_ = make_complex_from_top_faces(
-				maximal_faces_.begin(),
-				maximal_faces_.end(),
-				points_.begin(),
-				points_.end()
-		);
-	}
+  Complex& complex_;
+  typedef typename Complex::Vertex_handle Vertex_handle;
+  typedef typename Complex::Simplex_handle Simplex_handle;
+  typedef typename Complex::Point Point;
+
+  const bool load_only_points_;
+  std::vector<Point> points_;
+  std::vector<Simplex_handle> maximal_faces_;
+
+ public:
+  explicit Skeleton_blocker_off_visitor_reader(Complex& complex, bool load_only_points = false) :
+      complex_(complex),
+      load_only_points_(load_only_points) { }
+
+  void init(int dim, int num_vertices, int num_faces, int num_edges) {
+    maximal_faces_.reserve(num_faces);
+    points_.reserve(num_vertices);
+  }
+
+  void point(const std::vector<double>& point) {
+    points_.emplace_back(point.begin(), point.end());
+  }
+
+  void maximal_face(const std::vector<int>& face) {
+    if (!load_only_points_) {
+      Simplex_handle s;
+      for (auto x : face)
+        s.add_vertex(Vertex_handle(x));
+      maximal_faces_.emplace_back(s);
+    }
+  }
+
+  void done() {
+    complex_ = make_complex_from_top_faces(maximal_faces_.begin(), maximal_faces_.end(),
+                                           points_.begin(), points_.end() );
+  }
 };
 
-
 /**
  *@brief Class that allows to load a Skeleton_blocker_complex from an off file.
  */
 template<typename Complex>
 class Skeleton_blocker_off_reader {
-public:
-	/**
-	 * name_file : file to read
-	 * read_complex : complex that will receive the file content
-	 * read_only_points : specify true if only the points must be read
-	 */
-	Skeleton_blocker_off_reader(const std::string & name_file, Complex& read_complex, bool read_only_points = false,bool is_flag = false):valid_(false) {
-		std::ifstream stream(name_file);
-		if (stream.is_open()) {
-			if(is_flag || read_only_points){
-				Skeleton_blocker_off_flag_visitor_reader<Complex> off_visitor(read_complex, read_only_points);
-				Off_reader off_reader(stream);
-				valid_ = off_reader.read(off_visitor);
-			}
-			else{
-				Skeleton_blocker_off_visitor_reader<Complex> off_visitor(read_complex, read_only_points);
-				Off_reader off_reader(stream);
-				valid_ = off_reader.read(off_visitor);
-			}
-		}
-	}
-
-	/**
-	 * return true iff reading did not meet problems.
-	 */
-	bool is_valid() const {
-		return valid_;
-	}
-
-private:
-	bool valid_;
+ public:
+  /**
+   * name_file : file to read
+   * read_complex : complex that will receive the file content
+   * read_only_points : specify true if only the points must be read
+   */
+  Skeleton_blocker_off_reader(const std::string & name_file, Complex& read_complex,
+                              bool read_only_points = false, bool is_flag = false) : valid_(false) {
+    std::ifstream stream(name_file);
+    if (stream.is_open()) {
+      if (is_flag || read_only_points) {
+        Skeleton_blocker_off_flag_visitor_reader<Complex> off_visitor(read_complex, read_only_points);
+        Off_reader off_reader(stream);
+        valid_ = off_reader.read(off_visitor);
+      } else {
+        Skeleton_blocker_off_visitor_reader<Complex> off_visitor(read_complex, read_only_points);
+        Off_reader off_reader(stream);
+        valid_ = off_reader.read(off_visitor);
+      }
+    }
+  }
+
+  /**
+   * return true iff reading did not meet problems.
+   */
+  bool is_valid() const {
+    return valid_;
+  }
+
+ private:
+  bool valid_;
 };
 
-
 template<typename Complex>
-class Skeleton_blocker_off_writer{
-public:
-	/**
-	 * name_file : file where the off will be written
-	 * save_complex : complex that be outputted in the file
-	 * for now only save triangles.
-	 */
-	Skeleton_blocker_off_writer(const std::string & name_file, const Complex& save_complex){
-		typedef typename Complex::Vertex_handle Vertex_handle;
-
-		std::ofstream stream(name_file);
-		if (stream.is_open()) {
-			stream<<"OFF\n";
-			size_t num_triangles = std::distance(save_complex.triangle_range().begin(),save_complex.triangle_range().end());
-			stream<<save_complex.num_vertices()<<" "<<num_triangles<< " 0 \n";
-
-			//in case the complex has deactivated some vertices, eg only has vertices 0 2 5 7 for instance
-			//we compute a map from 0 2 5 7 to 0 1 2 3
-			std::map<Vertex_handle,size_t> vertex_num;
-			size_t current_vertex=0;
-
-			for(auto v : save_complex.vertex_range()){
-				vertex_num[v]=current_vertex++;
-				const auto& pt(save_complex.point(v));
-				for(auto x : pt)
-					stream<<x<<" ";
-				stream<<std::endl;
-			}
-
-			for(const auto & t : save_complex.triangle_range()){
-				stream<<"3 ";
-				for(auto x : t)
-					stream<<vertex_num[x]<<" ";
-				stream<<std::endl;
-			}
-			stream.close();
-		} else std::cerr<<"could not open file "<<name_file<<std::endl;
-	}
-
+class Skeleton_blocker_off_writer {
+ public:
+  /**
+   * name_file : file where the off will be written
+   * save_complex : complex that be outputted in the file
+   * for now only save triangles.
+   */
+  Skeleton_blocker_off_writer(const std::string & name_file, const Complex& save_complex) {
+    typedef typename Complex::Vertex_handle Vertex_handle;
+
+    std::ofstream stream(name_file);
+    if (stream.is_open()) {
+      stream << "OFF\n";
+      size_t num_triangles = std::distance(save_complex.triangle_range().begin(), save_complex.triangle_range().end());
+      stream << save_complex.num_vertices() << " " << num_triangles << " 0 \n";
+
+      // in case the complex has deactivated some vertices, eg only has vertices 0 2 5 7 for instance
+      // we compute a map from 0 2 5 7 to 0 1 2 3
+      std::map<Vertex_handle, size_t> vertex_num;
+      size_t current_vertex = 0;
+
+      for (auto v : save_complex.vertex_range()) {
+        vertex_num[v] = current_vertex++;
+        const auto& pt(save_complex.point(v));
+        for (auto x : pt)
+          stream << x << " ";
+        stream << std::endl;
+      }
+
+      for (const auto & t : save_complex.triangle_range()) {
+        stream << "3 ";
+        for (auto x : t)
+          stream << vertex_num[x] << " ";
+        stream << std::endl;
+      }
+      stream.close();
+    } else {
+      std::cerr << "could not open file " << name_file << std::endl;
+    }
+  }
 };
 
-
 }  // namespace skbl
 
-
 }  // namespace Gudhi
 
-
 #endif  // SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_SKELETON_BLOCKER_OFF_IO_H_
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
index 5b2c1f27..700830f2 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h
@@ -63,1479 +63,1466 @@ namespace skbl {
  */
 template<class SkeletonBlockerDS>
 class Skeleton_blocker_complex {
-	template<class ComplexType> friend class Vertex_iterator;
-	template<class ComplexType> friend class Neighbors_vertices_iterator;
-	template<class ComplexType> friend class Edge_iterator;
-	template<class ComplexType> friend class Edge_around_vertex_iterator;
-
-	template<class ComplexType> friend class Skeleton_blocker_link_complex;
-	template<class ComplexType> friend class Skeleton_blocker_link_superior;
-	template<class ComplexType> friend class Skeleton_blocker_sub_complex;
-
-public:
-	/**
-	 * @brief The type of stored vertex node, specified by the template SkeletonBlockerDS
-	 */
-	typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex;
-
-	/**
-	 * @brief The type of stored  edge node, specified by the template SkeletonBlockerDS
-	 */
-	typedef typename SkeletonBlockerDS::Graph_edge Graph_edge;
-
-	typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle;
-
-	/**
-	 * @brief The type of an handle to a vertex of the complex.
-	 */
-	typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle;
-	typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle;
-
-	/**
-	 * @brief A ordered set of integers that represents a simplex.
-	 */
-	typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle;
-	typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle;
-
-	/**
-	 * @brief Handle to a blocker of the complex.
-	 */
-	typedef Simplex_handle* Blocker_handle;
-
-	typedef typename Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_iterator;
-	typedef typename Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator;
-
-protected:
-	typedef typename boost::adjacency_list<boost::setS,  // edges
-			boost::vecS,  // vertices
-			boost::undirectedS, Graph_vertex, Graph_edge> Graph;
-	// todo/remark : edges are not sorted, it heavily penalizes computation for SuperiorLink
-	// (eg Link with greater vertices)
-	// that burdens simplex iteration / complex initialization via list of simplices.
-	// to avoid that, one should modify the graph by storing two lists of adjacency for every
-	// vertex, the one with superior and the one with lower vertices, that way there is
-	// no more extra cost for computation of SuperiorLink
-	typedef typename boost::graph_traits<Graph>::vertex_iterator boost_vertex_iterator;
-	typedef typename boost::graph_traits<Graph>::edge_iterator boost_edge_iterator;
-
-protected:
-	typedef typename boost::graph_traits<Graph>::adjacency_iterator boost_adjacency_iterator;
-
-public:
-	/**
-	 * @brief Handle to an edge of the complex.
-	 */
-	typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle;
-
-protected:
-	typedef std::multimap<Vertex_handle, Simplex_handle *> BlockerMap;
-	typedef typename std::multimap<Vertex_handle, Simplex_handle *>::value_type BlockerPair;
-	typedef typename std::multimap<Vertex_handle, Simplex_handle *>::iterator BlockerMapIterator;
-	typedef typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator BlockerMapConstIterator;
-
-protected:
-	int num_vertices_;
-	int num_blockers_;
-
-	typedef Skeleton_blocker_complex_visitor<Vertex_handle> Visitor;
-	// typedef Visitor* Visitor_ptr;
-	Visitor* visitor;
-
-	/**
-	 * @details If 'x' is a Vertex_handle of a vertex in the complex then degree[x] = d is its degree.
-	 *
-	 * This quantity is updated when adding/removing edge.
-	 *
-	 * This is useful because the operation
-	 * list.size() is done in linear time.
-	 */
-	std::vector<boost_vertex_handle> degree_;
-	Graph skeleton; /** 1-skeleton of the simplicial complex. */
-
-	/** Each vertex can access to the blockers passing through it. */
-	BlockerMap blocker_map_;
-
-public:
-	/////////////////////////////////////////////////////////////////////////////
-	/** @name Constructors, Destructors
-	 */
-	//@{
-
-	/**
-	 *@brief constructs a simplicial complex with a given number of vertices and a visitor.
-	 */
-	explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL)
-	: visitor(visitor_) {
-		clear();
-		for (int i = 0; i < num_vertices_; ++i) {
-			add_vertex();
-		}
-	}
-
-private:
-	//	typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie;
-	typedef Trie<Simplex_handle> STrie;
-
-public:
-	/**
-	 * @brief Constructor with a list of simplices.
-	 * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not).
-	 */
-	template<typename SimpleHandleOutputIterator>
-	Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, bool is_flag_complex = false, Visitor* visitor_ = NULL)
-	: num_vertices_(0),
-	  num_blockers_(0),
-	  visitor(visitor_) {
-		add_vertex_and_edges(simplex_begin,simplex_end);
-
-		if (!is_flag_complex)
-			// need to compute blockers
-			add_blockers(simplex_begin,simplex_end);
-	}
-
-private:
-	template<typename SimpleHandleOutputIterator>
-	void add_vertex_and_edges(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end){
-		std::vector<std::pair<Vertex_handle, Vertex_handle>> edges;
-		// first pass to add vertices and edges
-		int num_vertex = -1;
-		for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) {
-			if (s_it->dimension() == 0) num_vertex = (std::max)(num_vertex, s_it->first_vertex().vertex);
-			if (s_it->dimension() == 1) edges.emplace_back(s_it->first_vertex(), s_it->last_vertex());
-		}
-		while (num_vertex-- >= 0) add_vertex();
-
-		for (const auto& e : edges)
-			add_edge(e.first, e.second);
-	}
-
-	template<typename SimpleHandleOutputIterator>
-	void add_blockers(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end){
-		Tries<Simplex_handle> tries(num_vertices(),simplex_begin,simplex_end);
-		tries.init_next_dimension();
-		auto simplices(tries.next_dimension_simplices());
-
-		while(!simplices.empty()){
-			simplices = tries.next_dimension_simplices();
-			for(auto& sigma : simplices){
-				//common_positive_neighbors is the set of vertices u such that
-				//for all s in sigma, us is an edge and u>s
-				Simplex_handle common_positive_neighbors(tries.positive_neighbors(sigma.last_vertex()));
-				for(auto sigma_it = sigma.rbegin(); sigma_it != sigma.rend(); ++sigma_it/**/)
-					if(sigma_it != sigma.rbegin())
-						common_positive_neighbors.intersection(tries.positive_neighbors(*sigma_it));
-
-				for(auto x : common_positive_neighbors){
-					//first test that all edges sx are here for all s in sigma
-					bool all_edges_here = true;
-					for(auto s : sigma)
-						if(!contains_edge(x,s)){
-							all_edges_here = false;
-							break;
-						}
-					if(!all_edges_here) continue;
-
-					// all edges of sigma \cup x are here
-					//we have a blocker if all proper faces of sigma \cup x
-					// are in the complex and if sigma \cup x is not in the complex
-					//the first is equivalent at checking if blocks(sigma \cup x) is true
-					//as all blockers of lower dimension have already been computed
-					sigma.add_vertex(x);
-					if(!tries.contains(sigma)&&!blocks(sigma))
-						add_blocker(sigma);
-					sigma.remove_vertex(x);
-				}
-
-			}
-		}
-	}
-
-
-public:
-
-
-	// We cannot use the default copy constructor since we need
-	// to make a copy of each of the blockers
-
-	Skeleton_blocker_complex(const Skeleton_blocker_complex& copy) {
-		visitor = NULL;
-		degree_ = copy.degree_;
-		skeleton = Graph(copy.skeleton);
-		num_vertices_ = copy.num_vertices_;
-
-		num_blockers_ = 0;
-		// we copy the blockers
-		for (auto blocker : copy.const_blocker_range()) {
-			add_blocker(*blocker);
-		}
-	}
-
-	/**
-	 */
-	Skeleton_blocker_complex& operator=(const Skeleton_blocker_complex& copy) {
-		clear();
-		visitor = NULL;
-		degree_ = copy.degree_;
-		skeleton = Graph(copy.skeleton);
-		num_vertices_ = copy.num_vertices_;
-
-		num_blockers_ = 0;
-		// we copy the blockers
-		for (auto blocker : copy.const_blocker_range())
-			add_blocker(*blocker);
-		return *this;
-	}
-
-
-	/**
-	 * return true if both complexes have the same simplices.
-	 */
-	bool operator==(const Skeleton_blocker_complex& other) const{
-		if(other.num_vertices() != num_vertices()) return false;
-		if(other.num_edges() != num_edges()) return false;
-		if(other.num_blockers() != num_blockers()) return false;
-
-		for(auto v : vertex_range())
-			if(!other.contains_vertex(v)) return false;
-
-		for(auto e : edge_range())
-			if(!other.contains_edge(first_vertex(e),second_vertex(e))) return false;
-
-		for(const auto b : const_blocker_range())
-			if(!other.contains_blocker(*b)) return false;
-
-		return true;
-	}
-
-	bool operator!=(const Skeleton_blocker_complex& other) const{
-		return !(*this == other);
-	}
-
-	/**
-	 * The destructor delete all blockers allocated.
-	 */
-	virtual ~Skeleton_blocker_complex() {
-		clear();
-	}
-
-	/**
-	 * @details Clears the simplicial complex. After a call to this function,
-	 * blockers are destroyed. The 1-skeleton and the set of blockers
-	 * are both empty.
-	 */
-	virtual void clear() {
-		// xxx for now the responsabilty of freeing the visitor is for
-		// the user
-		visitor = NULL;
-
-		degree_.clear();
-		num_vertices_ = 0;
-
-		remove_blockers();
-
-		skeleton.clear();
-	}
-
-	/**
-	 *@brief allows to change the visitor.
-	 */
-	void set_visitor(Visitor* other_visitor) {
-		visitor = other_visitor;
-	}
-
-	//@}
-
-	/////////////////////////////////////////////////////////////////////////////
-	/** @name Vertices operations
-	 */
-	//@{
-public:
-	/**
-	 * @brief Return a local Vertex_handle of a vertex given a global one.
-	 * @remark Assume that the vertex is present in the complex.
-	 */
-	Vertex_handle operator[](Root_vertex_handle global) const {
-		auto local(get_address(global));
-		assert(local);
-		return *local;
-	}
-
-	/**
-	 * @brief Return the vertex node associated to local Vertex_handle.
-	 * @remark Assume that the vertex is present in the complex.
-	 */
-	Graph_vertex& operator[](Vertex_handle address) {
-		assert(
-				0 <= address.vertex && address.vertex < boost::num_vertices(skeleton));
-		return skeleton[address.vertex];
-	}
-
-	/**
-	 * @brief Return the vertex node associated to local Vertex_handle.
-	 * @remark Assume that the vertex is present in the complex.
-	 */
-	const Graph_vertex& operator[](Vertex_handle address) const {
-		assert(
-				0 <= address.vertex && address.vertex < boost::num_vertices(skeleton));
-		return skeleton[address.vertex];
-	}
-
-	/**
-	 * @brief Adds a vertex to the simplicial complex and returns its Vertex_handle.
-	 */
-	Vertex_handle add_vertex() {
-		Vertex_handle address(boost::add_vertex(skeleton));
-		num_vertices_++;
-		(*this)[address].activate();
-		// safe since we now that we are in the root complex and the field 'address' and 'id'
-		// are identical for every vertices
-		(*this)[address].set_id(Root_vertex_handle(address.vertex));
-		degree_.push_back(0);
-		if (visitor)
-			visitor->on_add_vertex(address);
-		return address;
-	}
-
-	/**
-	 * @brief Remove a vertex from the simplicial complex
-	 * @remark It just deactivates the vertex with a boolean flag but does not
-	 * remove it from vertices from complexity issues.
-	 */
-	void remove_vertex(Vertex_handle address) {
-		assert(contains_vertex(address));
-		// We remove b
-		boost::clear_vertex(address.vertex, skeleton);
-		(*this)[address].deactivate();
-		num_vertices_--;
-		degree_[address.vertex] = -1;
-		if (visitor)
-			visitor->on_remove_vertex(address);
-	}
-
-	/**
-	 */
-	bool contains_vertex(Vertex_handle u) const {
-		if (u.vertex < 0 || u.vertex >= boost::num_vertices(skeleton))
-			return false;
-		return (*this)[u].is_active();
-	}
-
-	/**
-	 */
-	bool contains_vertex(Root_vertex_handle u) const {
-		boost::optional<Vertex_handle> address = get_address(u);
-		return address && (*this)[*address].is_active();
-	}
-
-	/**
-	 * @return true iff the simplicial complex contains all vertices
-	 * of simplex sigma
-	 */
-	bool contains_vertices(const Simplex_handle & sigma) const {
-		for (auto vertex : sigma)
-			if (!contains_vertex(vertex))
-				return false;
-		return true;
-	}
-
-	/**
-	 * @brief Given an Id return the address of the vertex having this Id in the complex.
-	 * @remark For a simplicial complex, the address is the id but it may not be the case for a SubComplex.
-	 */
-	virtual boost::optional<Vertex_handle> get_address(
-			Root_vertex_handle id) const {
-		boost::optional<Vertex_handle> res;
-		if (id.vertex < boost::num_vertices(skeleton))
-			res = Vertex_handle(id.vertex);  // xxx
-		return res;
-	}
-
-	/**
-	 * return the id of a vertex of adress local present in the graph
-	 */
-	Root_vertex_handle get_id(Vertex_handle local) const {
-		assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton));
-		return (*this)[local].get_id();
-	}
-
-	/**
-	 * @brief Convert an address of a vertex of a complex to the address in
-	 * the current complex.
-	 * @details
-	 * If the current complex is a sub (or sup) complex of 'other', it converts
-	 * the address of a vertex v expressed in 'other' to the address of the vertex
-	 * v in the current one.
-	 * @remark this methods uses Root_vertex_handle to identify the vertex and
-	 * assumes the vertex is present in the current complex.
-	 */
-	Vertex_handle convert_handle_from_another_complex(
-			const Skeleton_blocker_complex& other, Vertex_handle vh_in_other) const {
-		auto vh_in_current_complex = get_address(other.get_id(vh_in_other));
-		assert(vh_in_current_complex);
-		return *vh_in_current_complex;
-	}
-
-	/**
-	 * @brief return the graph degree of a vertex.
-	 */
-	int degree(Vertex_handle local) const {
-		assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton));
-		return degree_[local.vertex];
-	}
-
-	//@}
-
-	/////////////////////////////////////////////////////////////////////////////
-	/** @name Edges operations
-	 */
-	//@{
-public:
-	/**
-	 * @brief return an edge handle if the two vertices forms
-	 * an edge in the complex
-	 */
-	boost::optional<Edge_handle> operator[](
-			const std::pair<Vertex_handle, Vertex_handle>& ab) const {
-		boost::optional<Edge_handle> res;
-		std::pair<Edge_handle, bool> edge_pair(
-				boost::edge(ab.first.vertex, ab.second.vertex, skeleton));
-		if (edge_pair.second)
-			res = edge_pair.first;
-		return res;
-	}
-
-	/**
-	 * @brief returns the stored node associated to an edge
-	 */
-	Graph_edge& operator[](Edge_handle edge_handle) {
-		return skeleton[edge_handle];
-	}
-
-	/**
-	 * @brief returns the stored node associated to an edge
-	 */
-	const Graph_edge& operator[](Edge_handle edge_handle) const {
-		return skeleton[edge_handle];
-	}
-
-	/**
-	 * @brief returns the first vertex of an edge
-	 * @details it assumes that the edge is present in the complex
-	 */
-	Vertex_handle first_vertex(Edge_handle edge_handle) const {
-		return static_cast<Vertex_handle> (source(edge_handle, skeleton));
-	}
-
-	/**
-	 * @brief returns the first vertex of an edge
-	 * @details it assumes that the edge is present in the complex
-	 */
-	Vertex_handle second_vertex(Edge_handle edge_handle) const {
-		return static_cast<Vertex_handle> (target(edge_handle, skeleton));
-	}
-
-	/**
-	 * @brief returns the simplex made with the two vertices of the edge
-	 * @details it assumes that the edge is present in the complex
-
-	 */
-	Simplex_handle get_vertices(Edge_handle edge_handle) const {
-		auto edge((*this)[edge_handle]);
-		return Simplex_handle((*this)[edge.first()], (*this)[edge.second()]);
-	}
-
-	/**
-	 * @brief Adds an edge between vertices a and b and all its cofaces.
-	 */
-	Edge_handle add_edge(Vertex_handle a, Vertex_handle b) {
-		assert(contains_vertex(a) && contains_vertex(b));
-		assert(a != b);
-
-		auto edge_handle((*this)[std::make_pair(a, b)]);
-		// std::pair<Edge_handle,bool> pair_descr_bool = (*this)[std::make_pair(a,b)];
-		// Edge_handle edge_descr;
-		// bool edge_present = pair_descr_bool.second;
-		if (!edge_handle) {
-			edge_handle = boost::add_edge(a.vertex, b.vertex, skeleton).first;
-			(*this)[*edge_handle].setId(get_id(a), get_id(b));
-			degree_[a.vertex]++;
-			degree_[b.vertex]++;
-			if (visitor)
-				visitor->on_add_edge(a, b);
-		}
-		return *edge_handle;
-	}
-
-	/**
-	 * @brief Adds all edges and their cofaces of a simplex to the simplicial complex.
-	 */
-	void add_edges(const Simplex_handle & sigma) {
-		Simplex_handle_iterator i, j;
-		for (i = sigma.begin(); i != sigma.end(); ++i)
-			for (j = i, j++; j != sigma.end(); ++j)
-				add_edge(*i, *j);
-	}
-
-	/**
-	 * @brief Removes an edge from the simplicial complex and all its cofaces.
-	 * @details returns the former Edge_handle representing the edge
-	 */
-	virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b) {
-		bool found;
-		Edge_handle edge;
-		tie(edge, found) = boost::edge(a.vertex, b.vertex, skeleton);
-		if (found) {
-			if (visitor)
-				visitor->on_remove_edge(a, b);
-			boost::remove_edge(a.vertex, b.vertex, skeleton);
-			degree_[a.vertex]--;
-			degree_[b.vertex]--;
-		}
-		return edge;
-	}
-
-	/**
-	 * @brief Removes edge and its cofaces from the simplicial complex.
-	 */
-	void remove_edge(Edge_handle edge) {
-		assert(contains_vertex(first_vertex(edge)));
-		assert(contains_vertex(second_vertex(edge)));
-		remove_edge(first_vertex(edge), second_vertex(edge));
-	}
-
-	/**
-	 * @brief The complex is reduced to its set of vertices.
-	 * All the edges and blockers are removed.
-	 */
-	void keep_only_vertices() {
-		remove_blockers();
-
-		for (auto u : vertex_range()) {
-			while (this->degree(u) > 0) {
-				Vertex_handle v(*(adjacent_vertices(u.vertex, this->skeleton).first));
-				this->remove_edge(u, v);
-			}
-		}
-	}
-
-	/**
-	 * @return true iff the simplicial complex contains an edge between
-	 * vertices a and b
-	 */
-	bool contains_edge(Vertex_handle a, Vertex_handle b) const {
-		// if (a.vertex<0 || b.vertex <0) return false;
-		return boost::edge(a.vertex, b.vertex, skeleton).second;
-	}
-
-	/**
-	 * @return true iff the simplicial complex contains all vertices
-	 * and all edges of simplex sigma
-	 */
-	bool contains_edges(const Simplex_handle & sigma) const {
-		for (auto i = sigma.begin(); i != sigma.end(); ++i) {
-			if (!contains_vertex(*i))
-				return false;
-			for (auto j = i; ++j != sigma.end();) {
-				if (!contains_edge(*i, *j))
-					return false;
-			}
-		}
-		return true;
-	}
-	//@}
-
-	/////////////////////////////////////////////////////////////////////////////
-	/** @name Blockers operations
-	 */
-	//@{
-
-	/**
-	 * @brief Adds the simplex to the set of blockers and
-	 * returns a Blocker_handle toward it if was not present before and 0 otherwise.
-	 */
-	Blocker_handle add_blocker(const Simplex_handle& blocker) {
-		assert(blocker.dimension() > 1);
-		if (contains_blocker(blocker)) {
-			// std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl;
-			return 0;
-		} else {
-			if (visitor)
-				visitor->on_add_blocker(blocker);
-			Blocker_handle blocker_pt = new Simplex_handle(blocker);
-			num_blockers_++;
-			auto vertex = blocker_pt->begin();
-			while (vertex != blocker_pt->end()) {
-				blocker_map_.insert(BlockerPair(*vertex, blocker_pt));
-				++vertex;
-			}
-			return blocker_pt;
-		}
-	}
-
-protected:
-	/**
-	 * @brief Adds the simplex to the set of blockers
-	 */
-	void add_blocker(Blocker_handle blocker) {
-		if (contains_blocker(*blocker)) {
-			// std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl;
-			return;
-		} else {
-			if (visitor)
-				visitor->on_add_blocker(*blocker);
-			num_blockers_++;
-			auto vertex = blocker->begin();
-			while (vertex != blocker->end()) {
-				blocker_map_.insert(BlockerPair(*vertex, blocker));
-				++vertex;
-			}
-		}
-	}
-
-protected:
-	/**
-	 * Removes sigma from the blocker map of vertex v
-	 */
-	void remove_blocker(const Blocker_handle sigma, Vertex_handle v) {
-		Complex_blocker_around_vertex_iterator blocker;
-		for (blocker = blocker_range(v).begin(); blocker != blocker_range(v).end();
-				++blocker) {
-			if (*blocker == sigma)
-				break;
-		}
-		if (*blocker != sigma) {
-			std::cerr
-			<< "bug ((*blocker).second == sigma) ie try to remove a blocker not present\n";
-			assert(false);
-		} else {
-			blocker_map_.erase(blocker.current_position());
-		}
-	}
-
-public:
-	/**
-	 * @brief Removes the simplex from the set of blockers.
-	 * @remark sigma has to belongs to the set of blockers
-	 */
-	void remove_blocker(const Blocker_handle sigma) {
-		for (auto vertex : *sigma)
-			remove_blocker(sigma, vertex);
-		num_blockers_--;
-	}
-
-	/**
-	 * @brief Remove all blockers, in other words, it expand the simplicial
-	 * complex to the smallest flag complex that contains it.
-	 */
-	void remove_blockers() {
-		// Desallocate the blockers
-		while (!blocker_map_.empty()) {
-			delete_blocker(blocker_map_.begin()->second);
-		}
-		num_blockers_ = 0;
-		blocker_map_.clear();
-	}
-
-protected:
-	/**
-	 * Removes the simplex sigma from the set of blockers.
-	 * sigma has to belongs to the set of blockers
-	 *
-	 * @remark contrarily to delete_blockers does not call the destructor
-	 */
-	void remove_blocker(const Simplex_handle& sigma) {
-		assert(contains_blocker(sigma));
-		for (auto vertex : sigma)
-			remove_blocker(sigma, vertex);
-		num_blockers_--;
-	}
-
-public:
-	/**
-	 * Removes the simplex s from the set of blockers
-	 * and desallocate s.
-	 */
-	void delete_blocker(Blocker_handle sigma) {
-		if (visitor)
-			visitor->on_delete_blocker(sigma);
-		remove_blocker(sigma);
-		delete sigma;
-	}
-
-	/**
-	 * @return true iff s is a blocker of the simplicial complex
-	 */
-	bool contains_blocker(const Blocker_handle s) const {
-		if (s->dimension() < 2)
-			return false;
-
-		Vertex_handle a = s->first_vertex();
-
-		for (const auto blocker : const_blocker_range(a)) {
-			if (s == *blocker)
-				return true;
-		}
-		return false;
-	}
-
-	/**
-	 * @return true iff s is a blocker of the simplicial complex
-	 */
-	bool contains_blocker(const Simplex_handle & s) const {
-		if (s.dimension() < 2)
-			return false;
-
-		Vertex_handle a = s.first_vertex();
-
-		for (auto blocker : const_blocker_range(a)) {
-			if (s == *blocker)
-				return true;
-		}
-		return false;
-	}
-
-private:
-	/**
-	 * @return true iff a blocker of the simplicial complex
-	 * is a face of sigma.
-	 */
-	bool blocks(const Simplex_handle & sigma) const {
-		for(auto s : sigma)
-			for (auto blocker : const_blocker_range(s))
-				if (sigma.contains(*blocker))
-					return true;
-		return false;
-
-	}
-
-	//@}
-
-protected:
-	/**
-	 * @details Adds to simplex the neighbours of v e.g. \f$ n \leftarrow n \cup N(v) \f$.
-	 * If keep_only_superior is true then only vertices that are greater than v are added.
-	 */
-	virtual void add_neighbours(Vertex_handle v, Simplex_handle & n,
-			bool keep_only_superior = false) const {
-		boost_adjacency_iterator ai, ai_end;
-		for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end;
-				++ai) {
-			if (keep_only_superior) {
-				if (*ai > v.vertex) {
-					n.add_vertex(Vertex_handle(*ai));
-				}
-			} else {
-				n.add_vertex(Vertex_handle(*ai));
-			}
-		}
-	}
-
-	/**
-	 * @details Add to simplex res all vertices which are
-	 * neighbours of alpha: ie \f$ res \leftarrow res \cup N(alpha) \f$.
-	 *
-	 * If 'keep_only_superior' is true then only vertices that are greater than alpha are added.
-	 * todo revoir
-	 *
-	 */
-	virtual void add_neighbours(const Simplex_handle &alpha, Simplex_handle & res,
-			bool keep_only_superior = false) const {
-		res.clear();
-		auto alpha_vertex = alpha.begin();
-		add_neighbours(*alpha_vertex, res, keep_only_superior);
-		for (alpha_vertex = (alpha.begin())++; alpha_vertex != alpha.end();
-				++alpha_vertex)
-			keep_neighbours(*alpha_vertex, res, keep_only_superior);
-	}
-
-	/**
-	 * @details Remove from simplex n all vertices which are
-	 * not neighbours of v e.g. \f$ res \leftarrow res \cap N(v) \f$.
-	 * If 'keep_only_superior' is true then only vertices that are greater than v are keeped.
-	 */
-	virtual void keep_neighbours(Vertex_handle v, Simplex_handle& res,
-			bool keep_only_superior = false) const {
-		Simplex_handle nv;
-		add_neighbours(v, nv, keep_only_superior);
-		res.intersection(nv);
-	}
-
-	/**
-	 * @details Remove from simplex all vertices which are
-	 * neighbours of v eg \f$ res \leftarrow res \setminus N(v) \f$.
-	 * If 'keep_only_superior' is true then only vertices that are greater than v are added.
-	 */
-	virtual void remove_neighbours(Vertex_handle v, Simplex_handle & res,
-			bool keep_only_superior = false) const {
-		Simplex_handle nv;
-		add_neighbours(v, nv, keep_only_superior);
-		res.difference(nv);
-	}
-
-public:
-	typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex> Link_complex;
-
-	/**
-	 * Constructs the link of 'simplex' with points coordinates.
-	 */
-	Link_complex link(Vertex_handle v) const {
-		return Link_complex(*this, Simplex_handle(v));
-	}
-
-	/**
-	 * Constructs the link of 'simplex' with points coordinates.
-	 */
-	Link_complex link(Edge_handle edge) const {
-		return Link_complex(*this, edge);
-	}
-
-	/**
-	 * Constructs the link of 'simplex' with points coordinates.
-	 */
-	Link_complex link(const Simplex_handle& simplex) const {
-		return Link_complex(*this, simplex);
-	}
-
-	/**
-	 * @brief Compute the local vertices of 's' in the current complex
-	 * If one of them is not present in the complex then the return value is uninitialized.
-	 *
-	 *
-	 */
-	// xxx rename get_address et place un using dans sub_complex
-
-	boost::optional<Simplex_handle> get_simplex_address(
-			const Root_simplex_handle& s) const {
-		boost::optional<Simplex_handle> res;
-
-		Simplex_handle s_address;
-		// Root_simplex_const_iterator i;
-		for (auto i = s.begin(); i != s.end(); ++i) {
-			boost::optional<Vertex_handle> address = get_address(*i);
-			if (!address)
-				return res;
-			else
-				s_address.add_vertex(*address);
-		}
-		res = s_address;
-		return res;
-	}
-
-	/**
-	 * @brief returns a simplex with vertices which are the id of vertices of the
-	 * argument.
-	 */
-	Root_simplex_handle get_id(const Simplex_handle& local_simplex) const {
-		Root_simplex_handle global_simplex;
-		for (auto x = local_simplex.begin(); x != local_simplex.end(); ++x) {
-			global_simplex.add_vertex(get_id(*x));
-		}
-		return global_simplex;
-	}
-
-	/**
-	 * @brief returns true iff the simplex s belongs to the simplicial
-	 * complex.
-	 */
-	virtual bool contains(const Simplex_handle & s) const {
-		if (s.dimension() == -1) {
-			return false;
-		} else if (s.dimension() == 0) {
-			return contains_vertex(s.first_vertex());
-		} else {
-			return (contains_edges(s) && !blocks(s));
-		}
-	}
-
-	/*
-	 * @brief returnrs true iff the complex is empty.
-	 */
-	bool empty() const {
-		return num_vertices() == 0;
-	}
-
-	/*
-	 * @brief returns the number of vertices in the complex.
-	 */
-	int num_vertices() const {
-		// remark boost::num_vertices(skeleton) counts deactivated vertices
-		return num_vertices_;
-	}
-
-	/*
-	 * @brief returns the number of edges in the complex.
-	 * @details currently in O(n)
-	 */
-	// todo cache the value
-
-	int num_edges() const {
-		return boost::num_edges(skeleton);
-	}
-
-	int num_triangles() const{
-		auto triangles = triangle_range();
-		return std::distance(triangles.begin(),triangles.end());
-	}
-	/*
-	 * @brief returns the number of blockers in the complex.
-	 */
-	int num_blockers() const {
-		return num_blockers_;
-	}
-
-	/*
-	 * @brief returns true iff the graph of the 1-skeleton of the complex is complete.
-	 */
-	bool complete() const {
-		return (num_vertices() * (num_vertices() - 1)) / 2 == num_edges();
-	}
-
-	/**
-	 * @brief returns the number of connected components in the graph of the 1-skeleton.
-	 */
-	int num_connected_components() const {
-		int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices();
-		std::vector<int> component(skeleton.vertex_set().size());
-		return boost::connected_components(this->skeleton, &component[0])
-		- num_vert_collapsed;
-	}
-
-	/**
-	 * @brief %Test if the complex is a cone.
-	 * @details Runs in O(n) where n is the number of vertices.
-	 */
-	bool is_cone() const {
-		if (num_vertices() == 0)
-			return false;
-		if (num_vertices() == 1)
-			return true;
-		for (auto vi : vertex_range()) {
-			// xxx todo faire une methode bool is_in_blocker(Vertex_handle)
-			if (blocker_map_.find(vi) == blocker_map_.end()) {
-				// no blocker passes through the vertex, we just need to
-				// check if the current vertex is linked to all others vertices of the complex
-				if (degree_[vi.vertex] == num_vertices() - 1)
-					return true;
-			}
-		}
-		return false;
-	}
-
-	//@}
-	/** @Simplification operations
-	 */
-	//@{
-
-	/**
-	 * Returns true iff the blocker 'sigma' is popable.
-	 * To define popable, let us call 'L' the complex that
-	 * consists in the current complex without the blocker 'sigma'.
-	 * A blocker 'sigma' is then "popable" if the link of 'sigma'
-	 * in L is reducible.
-	 *
-	 */
-	bool is_popable_blocker(Blocker_handle sigma) const;
-
-	/**
-	 * Removes all the popable blockers of the complex and delete them.
-	 * @returns the number of popable blockers deleted
-	 */
-	void remove_popable_blockers();
-
-	/**
-	 * Removes all the popable blockers of the complex passing through v and delete them.
-	 */
-	void remove_popable_blockers(Vertex_handle v);
-
-	/**
-	 * @brief Removes all the popable blockers of the complex passing through v and delete them.
-	 * Also remove popable blockers in the neighborhood if they became popable.
-	 *
-	 */
-	void remove_all_popable_blockers(Vertex_handle v);
-
-	/**
-	 * Remove the star of the vertex 'v'
-	 */
-	void remove_star(Vertex_handle v) ;
-
-private:
-	/**
-	 * after removing the star of a simplex, blockers sigma that contains this simplex must be removed.
-	 * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added
-	 * whenever the dimension of tau is at least 2.
-	 */
-	void update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed);
-
-public:
-	/**
-	 * Remove the star of the edge connecting vertices a and b.
-	 * @returns the number of blocker that have been removed
-	 */
-	void remove_star(Vertex_handle a, Vertex_handle b) ;
-
-	/**
-	 * Remove the star of the edge 'e'.
-	 */
-	void remove_star(Edge_handle e) ;
-
-	/**
-	 * Remove the star of the simplex 'sigma' which needs to belong to the complex
-	 */
-	void remove_star(const Simplex_handle& sigma);
-
-	/**
-	 * @brief add a maximal simplex plus all its cofaces.
-	 * @details the simplex must have dimension greater than one (otherwise use add_vertex or add_edge).
-	 */
-	void add_simplex(const Simplex_handle& sigma);
-
-private:
-	/**
-	 * remove all blockers that contains sigma
-	 */
-	void remove_blocker_containing_simplex(const Simplex_handle& sigma) ;
-
-	/**
-	 * remove all blockers that contains sigma
-	 */
-	void remove_blocker_include_in_simplex(const Simplex_handle& sigma);
-
-public:
-	enum simplifiable_status {
-		NOT_HOMOTOPY_EQ, MAYBE_HOMOTOPY_EQ, HOMOTOPY_EQ
-	};
-
-	simplifiable_status is_remove_star_homotopy_preserving(const Simplex_handle& simplex) {
-		// todo write a virtual method 'link' in Skeleton_blocker_complex which will be overloaded by the current one of Skeleton_blocker_geometric_complex
-		// then call it there to build the link and return the value of link.is_contractible()
-		return MAYBE_HOMOTOPY_EQ;
-	}
-
-	enum contractible_status {
-		NOT_CONTRACTIBLE, MAYBE_CONTRACTIBLE, CONTRACTIBLE
-	};
-
-	/**
-	 * @brief %Test if the complex is reducible using a strategy defined in the class
-	 * (by default it tests if the complex is a cone)
-	 * @details Note that NO could be returned if some invariant ensures that the complex
-	 * is not a point (for instance if the euler characteristic is different from 1).
-	 * This function will surely have to return MAYBE in some case because the
-	 * associated problem is undecidable but it in practice, it can often
-	 * be solved with the help of geometry.
-	 */
-	virtual contractible_status is_contractible() const {
-		if (this->is_cone()) {
-			return CONTRACTIBLE;
-		} else {
-			return MAYBE_CONTRACTIBLE;
-		}
-	}
-	//@}
-
-	/** @Edge contraction operations
-	 */
-	//@{
-
-	/**
-	 * @return If ignore_popable_blockers is true
-	 * then the result is true iff the link condition at edge ab is satisfied
-	 * or equivalently iff no blocker contains ab.
-	 * If ignore_popable_blockers is false then the
-	 * result is true iff all blocker containing ab are popable.
-	 */
-	bool link_condition(Vertex_handle a, Vertex_handle b, bool ignore_popable_blockers = false) const{
-		for (auto blocker : this->const_blocker_range(a))
-			if (blocker->contains(b)) {
-				// false if ignore_popable_blockers is false
-				// otherwise the blocker has to be popable
-				return ignore_popable_blockers && is_popable_blocker(blocker);
-			}
-		return true;
-
-	}
-
-	/**
-	 * @return If ignore_popable_blockers is true
-	 * then the result is true iff the link condition at edge ab is satisfied
-	 * or equivalently iff no blocker contains ab.
-	 * If ignore_popable_blockers is false then the
-	 * result is true iff all blocker containing ab are popable.
-	 */
-	bool link_condition(Edge_handle e, bool ignore_popable_blockers = false) const {
-		const Graph_edge& edge = (*this)[e];
-		assert(this->get_address(edge.first()));
-		assert(this->get_address(edge.second()));
-		Vertex_handle a(*this->get_address(edge.first()));
-		Vertex_handle b(*this->get_address(edge.second()));
-		return link_condition(a, b, ignore_popable_blockers);
-	}
-
-protected:
-	/**
-	 * Compute simplices beta such that a.beta is an order 0 blocker
-	 * that may be used to construct a new blocker after contracting ab.
-	 * It requires that the link condition is satisfied.
-	 */
-	void tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const;
-
-private:
-	/**
-	 * @brief "Replace" the edge 'bx' by the edge 'ax'.
-	 * Assume that the edge 'bx' was present whereas 'ax' was not.
-	 * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'.
-	 * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added
-	 * and just before edge 'bx' had been removed. That way, it can
-	 * eventually access to information of 'ax'.
-	 */
-	void swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x);
-
-private:
-	/**
-	 * @brief removes all blockers passing through the edge 'ab'
-	 */
-	void delete_blockers_around_vertex(Vertex_handle v);
-
-	/**
-	 * @brief removes all blockers passing through the edge 'ab'
-	 */
-	void delete_blockers_around_edge(Vertex_handle a, Vertex_handle b);
-
-public:
-	/**
-	 * Contracts the edge.
-	 * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
-	 * it removes first all blockers passing through 'ab'
-	 */
-	void contract_edge(Edge_handle edge) {
-		contract_edge(this->first_vertex(edge), this->second_vertex(edge));
-	}
-
-	/**
-	 * Contracts the edge connecting vertices a and b.
-	 * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
-	 * it removes first all blockers passing through 'ab'
-	 */
-	void contract_edge(Vertex_handle a, Vertex_handle b);
-
-private:
-	void get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, std::set<Simplex_handle>& blockers_to_add);
-
-	/**
-	 * delete all blockers that passes through a or b
-	 */
-	void delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b);
-	void update_edges_after_contraction(Vertex_handle a, Vertex_handle b) ;
-
-	void notify_changed_edges(Vertex_handle a) ;
-	//@}
-
-public:
-	/////////////////////////////////////////////////////////////////////////////
-	/** @name Vertex iterators
-	 */
-	//@{
-	typedef Vertex_iterator<Skeleton_blocker_complex> Complex_vertex_iterator;
-
-	//
-	// Range over the vertices of the simplicial complex.
-	// Methods .begin() and .end() return a Complex_vertex_iterator.
-	//
-	typedef boost::iterator_range<Complex_vertex_iterator> Complex_vertex_range;
-
-	/**
-	 * @brief Returns a Complex_vertex_range over all vertices of the complex
-	 */
-	Complex_vertex_range vertex_range() const {
-		auto begin = Complex_vertex_iterator(this);
-		auto end = Complex_vertex_iterator(this, 0);
-		return Complex_vertex_range(begin, end);
-	}
-
-	typedef Neighbors_vertices_iterator<Skeleton_blocker_complex> Complex_neighbors_vertices_iterator;
-
-
-	typedef boost::iterator_range<Complex_neighbors_vertices_iterator> Complex_neighbors_vertices_range;
-
-	/**
-	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v
-	 */
-	Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const {
-		auto begin = Complex_neighbors_vertices_iterator(this, v);
-		auto end = Complex_neighbors_vertices_iterator(this, v, 0);
-		return Complex_neighbors_vertices_range(begin, end);
-	}
-
-	//@}
-
-	/** @name Edge iterators
-	 */
-	//@{
-
-	typedef Edge_iterator<Skeleton_blocker_complex> Complex_edge_iterator;
-
-
-	typedef boost::iterator_range<Complex_edge_iterator> Complex_edge_range;
-
-	/**
-	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex
-	 */
-	Complex_edge_range edge_range() const {
-		auto begin = Complex_edge_iterator(this);
-		auto end = Complex_edge_iterator(this, 0);
-		return Complex_edge_range(begin, end);
-	}
-
-
-	typedef Edge_around_vertex_iterator<Skeleton_blocker_complex> Complex_edge_around_vertex_iterator;
-
-
-	typedef boost::iterator_range <Complex_edge_around_vertex_iterator>	Complex_edge_around_vertex_range;
-
-	/**
-	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes
-	 * through 'v'
-	 */
-	Complex_edge_around_vertex_range edge_range(Vertex_handle v) const {
-		auto begin = Complex_edge_around_vertex_iterator(this, v);
-		auto end = Complex_edge_around_vertex_iterator(this, v, 0);
-		return Complex_edge_around_vertex_range(begin, end);
-	}
-
-	//@}
-
-	/** @name Triangles iterators
-	 */
-	//@{
-private:
-	typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link;
-	typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link;
-
-public:
-	typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex, Superior_link> Superior_triangle_around_vertex_iterator;
-
-	typedef boost::iterator_range < Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link> > Complex_triangle_around_vertex_range;
-
-	/**
-	 * @brief Range over triangles around a vertex of the simplicial complex.
-	 * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
-	 *
-	 */
-	Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const {
-		auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v);
-		auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v, 0);
-		return Complex_triangle_around_vertex_range(begin, end);
-	}
-
-	typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range;
-
-
-	typedef Triangle_iterator<Skeleton_blocker_complex> Complex_triangle_iterator;
-
-
-	/**
-	 * @brief Range over triangles of the simplicial complex.
-	 * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
-	 *
-	 */
-	Complex_triangle_range triangle_range() const {
-		auto end = Triangle_iterator<Skeleton_blocker_complex>(this, 0);
-		if (empty()) {
-			return Complex_triangle_range(end, end);
-		} else {
-			auto begin = Triangle_iterator<Skeleton_blocker_complex>(this);
-			return Complex_triangle_range(begin, end);
-		}
-	}
-
-	//@}
-
-	/** @name Simplices iterators
-	 */
-	//@{
-	typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex, Link> Complex_simplex_around_vertex_iterator;
-
-	/**
-	 * @brief Range over the simplices of the simplicial complex around a vertex.
-	 * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator.
-	 */
-	typedef boost::iterator_range < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range;
-
-	/**
-	 * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex
-	 */
-	Complex_simplex_around_vertex_range simplex_range(Vertex_handle v) const {
-		assert(contains_vertex(v));
-		return Complex_simplex_around_vertex_range(
-				Complex_simplex_around_vertex_iterator(this, v),
-				Complex_simplex_around_vertex_iterator(this, v, true));
-	}
-
-	// typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator;
-	typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator;
-
-	typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range;
-
-	/**
-	 * @brief Returns a Complex_simplex_range over all the simplices of the complex
-	 */
-	Complex_simplex_range simplex_range() const {
-		Complex_simplex_iterator end(this, true);
-		if (empty()) {
-			return Complex_simplex_range(end, end);
-		} else {
-			Complex_simplex_iterator begin(this);
-			return Complex_simplex_range(begin, end);
-		}
-	}
-
-	//@}
-
-	/** @name Blockers iterators
-	 */
-	//@{
-private:
-	/**
-	 * @brief Iterator over the blockers adjacent to a vertex
-	 */
-	typedef Blocker_iterator_around_vertex_internal<
-			typename std::multimap<Vertex_handle, Simplex_handle *>::iterator,
-			Blocker_handle>
-	Complex_blocker_around_vertex_iterator;
-
-	/**
-	 * @brief Iterator over (constant) blockers adjacent to a vertex
-	 */
-	typedef Blocker_iterator_around_vertex_internal<
-			typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator,
-			const Blocker_handle>
-	Const_complex_blocker_around_vertex_iterator;
-
-	typedef boost::iterator_range <Complex_blocker_around_vertex_iterator> Complex_blocker_around_vertex_range;
-	typedef boost::iterator_range <Const_complex_blocker_around_vertex_iterator> Const_complex_blocker_around_vertex_range;
-
-public:
-	/**
-	 * @brief Returns a range of the blockers of the complex passing through a vertex
-	 */
-	Complex_blocker_around_vertex_range blocker_range(Vertex_handle v) {
-		auto begin = Complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v));
-		auto end = Complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v));
-		return Complex_blocker_around_vertex_range(begin, end);
-	}
-
-	/**
-	 * @brief Returns a range of the blockers of the complex passing through a vertex
-	 */
-	Const_complex_blocker_around_vertex_range const_blocker_range(Vertex_handle v) const {
-		auto begin = Const_complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v));
-		auto end = Const_complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v));
-		return Const_complex_blocker_around_vertex_range(begin, end);
-	}
-
-private:
-	/**
-	 * @brief Iterator over the blockers.
-	 */
-	typedef Blocker_iterator_internal<
-			typename std::multimap<Vertex_handle, Simplex_handle *>::iterator,
-			Blocker_handle>
-	Complex_blocker_iterator;
-
-	/**
-	 * @brief Iterator over the (constant) blockers.
-	 */
-	typedef Blocker_iterator_internal<
-			typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator,
-			const Blocker_handle>
-	Const_complex_blocker_iterator;
-
-	typedef boost::iterator_range <Complex_blocker_iterator> Complex_blocker_range;
-	typedef boost::iterator_range <Const_complex_blocker_iterator> Const_complex_blocker_range;
-
-public:
-	/**
-	 * @brief Returns a range of the blockers of the complex
-	 */
-	Complex_blocker_range blocker_range() {
-		auto begin = Complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end());
-		auto end = Complex_blocker_iterator(blocker_map_.end(), blocker_map_.end());
-		return Complex_blocker_range(begin, end);
-	}
-
-	/**
-	 * @brief Returns a range of the blockers of the complex
-	 */
-	Const_complex_blocker_range const_blocker_range() const {
-		auto begin = Const_complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end());
-		auto end = Const_complex_blocker_iterator(blocker_map_.end(), blocker_map_.end());
-		return Const_complex_blocker_range(begin, end);
-	}
-
-	//@}
-
-	/////////////////////////////////////////////////////////////////////////////
-	/** @name Print and IO methods
-	 */
-	//@{
-public:
-	std::string to_string() const {
-		std::ostringstream stream;
-		stream << num_vertices() << " vertices:\n" << vertices_to_string() << std::endl;
-		stream << num_edges() << " edges:\n" << edges_to_string() << std::endl;
-		stream << num_blockers() << " blockers:\n" << blockers_to_string() << std::endl;
-		return stream.str();
-	}
-
-	std::string vertices_to_string() const {
-		std::ostringstream stream;
-		for (auto vertex : vertex_range()) {
-			stream << "{" << (*this)[vertex].get_id() << "} ";
-		}
-		stream << std::endl;
-		return stream.str();
-	}
-
-	std::string edges_to_string() const {
-		std::ostringstream stream;
-		for (auto edge : edge_range())
-			stream << "{" << (*this)[edge].first() << "," << (*this)[edge].second() << "} ";
-		stream << std::endl;
-		return stream.str();
-	}
-
-	std::string blockers_to_string() const {
-		std::ostringstream stream;
-
-		for (auto b : const_blocker_range())
-			stream << *b << std::endl;
-		return stream.str();
-	}
-	//@}
-
-
-
-
-
+  template<class ComplexType> friend class Vertex_iterator;
+  template<class ComplexType> friend class Neighbors_vertices_iterator;
+  template<class ComplexType> friend class Edge_iterator;
+  template<class ComplexType> friend class Edge_around_vertex_iterator;
+
+  template<class ComplexType> friend class Skeleton_blocker_link_complex;
+  template<class ComplexType> friend class Skeleton_blocker_link_superior;
+  template<class ComplexType> friend class Skeleton_blocker_sub_complex;
+
+ public:
+  /**
+   * @brief The type of stored vertex node, specified by the template SkeletonBlockerDS
+   */
+  typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex;
+
+  /**
+   * @brief The type of stored  edge node, specified by the template SkeletonBlockerDS
+   */
+  typedef typename SkeletonBlockerDS::Graph_edge Graph_edge;
+
+  typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle;
+
+  /**
+   * @brief The type of an handle to a vertex of the complex.
+   */
+  typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle;
+  typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle;
+
+  /**
+   * @brief A ordered set of integers that represents a simplex.
+   */
+  typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle;
+  typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle;
+
+  /**
+   * @brief Handle to a blocker of the complex.
+   */
+  typedef Simplex_handle* Blocker_handle;
+
+  typedef typename Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_iterator;
+  typedef typename Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator;
+
+ protected:
+  typedef typename boost::adjacency_list<boost::setS,  // edges
+  boost::vecS,  // vertices
+  boost::undirectedS, Graph_vertex, Graph_edge> Graph;
+  // todo/remark : edges are not sorted, it heavily penalizes computation for SuperiorLink
+  // (eg Link with greater vertices)
+  // that burdens simplex iteration / complex initialization via list of simplices.
+  // to avoid that, one should modify the graph by storing two lists of adjacency for every
+  // vertex, the one with superior and the one with lower vertices, that way there is
+  // no more extra cost for computation of SuperiorLink
+  typedef typename boost::graph_traits<Graph>::vertex_iterator boost_vertex_iterator;
+  typedef typename boost::graph_traits<Graph>::edge_iterator boost_edge_iterator;
+
+ protected:
+  typedef typename boost::graph_traits<Graph>::adjacency_iterator boost_adjacency_iterator;
+
+ public:
+  /**
+   * @brief Handle to an edge of the complex.
+   */
+  typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle;
+
+ protected:
+  typedef std::multimap<Vertex_handle, Simplex_handle *> BlockerMap;
+  typedef typename std::multimap<Vertex_handle, Simplex_handle *>::value_type BlockerPair;
+  typedef typename std::multimap<Vertex_handle, Simplex_handle *>::iterator BlockerMapIterator;
+  typedef typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator BlockerMapConstIterator;
+
+ protected:
+  int num_vertices_;
+  int num_blockers_;
+
+  typedef Skeleton_blocker_complex_visitor<Vertex_handle> Visitor;
+  // typedef Visitor* Visitor_ptr;
+  Visitor* visitor;
+
+  /**
+   * @details If 'x' is a Vertex_handle of a vertex in the complex then degree[x] = d is its degree.
+   *
+   * This quantity is updated when adding/removing edge.
+   *
+   * This is useful because the operation
+   * list.size() is done in linear time.
+   */
+  std::vector<boost_vertex_handle> degree_;
+  Graph skeleton; /** 1-skeleton of the simplicial complex. */
+
+  /** Each vertex can access to the blockers passing through it. */
+  BlockerMap blocker_map_;
+
+ public:
+  /////////////////////////////////////////////////////////////////////////////
+  /** @name Constructors, Destructors
+   */
+  //@{
+
+  /**
+   *@brief constructs a simplicial complex with a given number of vertices and a visitor.
+   */
+  explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL)
+      : visitor(visitor_) {
+    clear();
+    for (int i = 0; i < num_vertices_; ++i) {
+      add_vertex();
+    }
+  }
+
+ private:
+  // typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie;
+  typedef Trie<Simplex_handle> STrie;
+
+ public:
+  /**
+   * @brief Constructor with a list of simplices.
+   * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not).
+   */
+  template<typename SimpleHandleOutputIterator>
+  Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end,
+                           bool is_flag_complex = false, Visitor* visitor_ = NULL)
+      : num_vertices_(0),
+      num_blockers_(0),
+      visitor(visitor_) {
+    add_vertex_and_edges(simplex_begin, simplex_end);
+
+    if (!is_flag_complex)
+      // need to compute blockers
+      add_blockers(simplex_begin, simplex_end);
+  }
+
+ private:
+  template<typename SimpleHandleOutputIterator>
+  void add_vertex_and_edges(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) {
+    std::vector<std::pair<Vertex_handle, Vertex_handle>> edges;
+    // first pass to add vertices and edges
+    int num_vertex = -1;
+    for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) {
+      if (s_it->dimension() == 0) num_vertex = (std::max)(num_vertex, s_it->first_vertex().vertex);
+      if (s_it->dimension() == 1) edges.emplace_back(s_it->first_vertex(), s_it->last_vertex());
+    }
+    while (num_vertex-- >= 0) add_vertex();
+
+    for (const auto& e : edges)
+      add_edge(e.first, e.second);
+  }
+
+  template<typename SimpleHandleOutputIterator>
+  void add_blockers(SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) {
+    Tries<Simplex_handle> tries(num_vertices(), simplex_begin, simplex_end);
+    tries.init_next_dimension();
+    auto simplices(tries.next_dimension_simplices());
+
+    while (!simplices.empty()) {
+      simplices = tries.next_dimension_simplices();
+      for (auto& sigma : simplices) {
+        // common_positive_neighbors is the set of vertices u such that
+        // for all s in sigma, us is an edge and u>s
+        Simplex_handle common_positive_neighbors(tries.positive_neighbors(sigma.last_vertex()));
+        for (auto sigma_it = sigma.rbegin(); sigma_it != sigma.rend(); ++sigma_it)
+          if (sigma_it != sigma.rbegin())
+            common_positive_neighbors.intersection(tries.positive_neighbors(*sigma_it));
+
+        for (auto x : common_positive_neighbors) {
+          // first test that all edges sx are here for all s in sigma
+          bool all_edges_here = true;
+          for (auto s : sigma)
+            if (!contains_edge(x, s)) {
+              all_edges_here = false;
+              break;
+            }
+          if (!all_edges_here) continue;
+
+          // all edges of sigma \cup x are here
+          // we have a blocker if all proper faces of sigma \cup x
+          // are in the complex and if sigma \cup x is not in the complex
+          // the first is equivalent at checking if blocks(sigma \cup x) is true
+          // as all blockers of lower dimension have already been computed
+          sigma.add_vertex(x);
+          if (!tries.contains(sigma) && !blocks(sigma))
+            add_blocker(sigma);
+          sigma.remove_vertex(x);
+        }
+      }
+    }
+  }
+
+ public:
+  // We cannot use the default copy constructor since we need
+  // to make a copy of each of the blockers
+
+  Skeleton_blocker_complex(const Skeleton_blocker_complex& copy) {
+    visitor = NULL;
+    degree_ = copy.degree_;
+    skeleton = Graph(copy.skeleton);
+    num_vertices_ = copy.num_vertices_;
+
+    num_blockers_ = 0;
+    // we copy the blockers
+    for (auto blocker : copy.const_blocker_range()) {
+      add_blocker(*blocker);
+    }
+  }
+
+  /**
+   */
+  Skeleton_blocker_complex& operator=(const Skeleton_blocker_complex& copy) {
+    clear();
+    visitor = NULL;
+    degree_ = copy.degree_;
+    skeleton = Graph(copy.skeleton);
+    num_vertices_ = copy.num_vertices_;
+
+    num_blockers_ = 0;
+    // we copy the blockers
+    for (auto blocker : copy.const_blocker_range())
+      add_blocker(*blocker);
+    return *this;
+  }
+
+  /**
+   * return true if both complexes have the same simplices.
+   */
+  bool operator==(const Skeleton_blocker_complex& other) const {
+    if (other.num_vertices() != num_vertices()) return false;
+    if (other.num_edges() != num_edges()) return false;
+    if (other.num_blockers() != num_blockers()) return false;
+
+    for (auto v : vertex_range())
+      if (!other.contains_vertex(v)) return false;
+
+    for (auto e : edge_range())
+      if (!other.contains_edge(first_vertex(e), second_vertex(e))) return false;
+
+    for (const auto b : const_blocker_range())
+      if (!other.contains_blocker(*b)) return false;
+
+    return true;
+  }
+
+  bool operator!=(const Skeleton_blocker_complex& other) const {
+    return !(*this == other);
+  }
+
+  /**
+   * The destructor delete all blockers allocated.
+   */
+  virtual ~Skeleton_blocker_complex() {
+    clear();
+  }
+
+  /**
+   * @details Clears the simplicial complex. After a call to this function,
+   * blockers are destroyed. The 1-skeleton and the set of blockers
+   * are both empty.
+   */
+  virtual void clear() {
+    // xxx for now the responsabilty of freeing the visitor is for
+    // the user
+    visitor = NULL;
+
+    degree_.clear();
+    num_vertices_ = 0;
+
+    remove_blockers();
+
+    skeleton.clear();
+  }
+
+  /**
+   *@brief allows to change the visitor.
+   */
+  void set_visitor(Visitor* other_visitor) {
+    visitor = other_visitor;
+  }
+
+  //@}
+
+  /////////////////////////////////////////////////////////////////////////////
+  /** @name Vertices operations
+   */
+  //@{
+ public:
+  /**
+   * @brief Return a local Vertex_handle of a vertex given a global one.
+   * @remark Assume that the vertex is present in the complex.
+   */
+  Vertex_handle operator[](Root_vertex_handle global) const {
+    auto local(get_address(global));
+    assert(local);
+    return *local;
+  }
+
+  /**
+   * @brief Return the vertex node associated to local Vertex_handle.
+   * @remark Assume that the vertex is present in the complex.
+   */
+  Graph_vertex& operator[](Vertex_handle address) {
+    assert(
+           0 <= address.vertex && address.vertex < boost::num_vertices(skeleton));
+    return skeleton[address.vertex];
+  }
+
+  /**
+   * @brief Return the vertex node associated to local Vertex_handle.
+   * @remark Assume that the vertex is present in the complex.
+   */
+  const Graph_vertex& operator[](Vertex_handle address) const {
+    assert(
+           0 <= address.vertex && address.vertex < boost::num_vertices(skeleton));
+    return skeleton[address.vertex];
+  }
+
+  /**
+   * @brief Adds a vertex to the simplicial complex and returns its Vertex_handle.
+   */
+  Vertex_handle add_vertex() {
+    Vertex_handle address(boost::add_vertex(skeleton));
+    num_vertices_++;
+    (*this)[address].activate();
+    // safe since we now that we are in the root complex and the field 'address' and 'id'
+    // are identical for every vertices
+    (*this)[address].set_id(Root_vertex_handle(address.vertex));
+    degree_.push_back(0);
+    if (visitor)
+      visitor->on_add_vertex(address);
+    return address;
+  }
+
+  /**
+   * @brief Remove a vertex from the simplicial complex
+   * @remark It just deactivates the vertex with a boolean flag but does not
+   * remove it from vertices from complexity issues.
+   */
+  void remove_vertex(Vertex_handle address) {
+    assert(contains_vertex(address));
+    // We remove b
+    boost::clear_vertex(address.vertex, skeleton);
+    (*this)[address].deactivate();
+    num_vertices_--;
+    degree_[address.vertex] = -1;
+    if (visitor)
+      visitor->on_remove_vertex(address);
+  }
+
+  /**
+   */
+  bool contains_vertex(Vertex_handle u) const {
+    if (u.vertex < 0 || u.vertex >= boost::num_vertices(skeleton))
+      return false;
+    return (*this)[u].is_active();
+  }
+
+  /**
+   */
+  bool contains_vertex(Root_vertex_handle u) const {
+    boost::optional<Vertex_handle> address = get_address(u);
+    return address && (*this)[*address].is_active();
+  }
+
+  /**
+   * @return true iff the simplicial complex contains all vertices
+   * of simplex sigma
+   */
+  bool contains_vertices(const Simplex_handle & sigma) const {
+    for (auto vertex : sigma)
+      if (!contains_vertex(vertex))
+        return false;
+    return true;
+  }
+
+  /**
+   * @brief Given an Id return the address of the vertex having this Id in the complex.
+   * @remark For a simplicial complex, the address is the id but it may not be the case for a SubComplex.
+   */
+  virtual boost::optional<Vertex_handle> get_address(
+                                                     Root_vertex_handle id) const {
+    boost::optional<Vertex_handle> res;
+    if (id.vertex < boost::num_vertices(skeleton))
+      res = Vertex_handle(id.vertex);  // xxx
+    return res;
+  }
+
+  /**
+   * return the id of a vertex of adress local present in the graph
+   */
+  Root_vertex_handle get_id(Vertex_handle local) const {
+    assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton));
+    return (*this)[local].get_id();
+  }
+
+  /**
+   * @brief Convert an address of a vertex of a complex to the address in
+   * the current complex.
+   * @details
+   * If the current complex is a sub (or sup) complex of 'other', it converts
+   * the address of a vertex v expressed in 'other' to the address of the vertex
+   * v in the current one.
+   * @remark this methods uses Root_vertex_handle to identify the vertex and
+   * assumes the vertex is present in the current complex.
+   */
+  Vertex_handle convert_handle_from_another_complex(const Skeleton_blocker_complex& other,
+                                                    Vertex_handle vh_in_other) const {
+    auto vh_in_current_complex = get_address(other.get_id(vh_in_other));
+    assert(vh_in_current_complex);
+    return *vh_in_current_complex;
+  }
+
+  /**
+   * @brief return the graph degree of a vertex.
+   */
+  int degree(Vertex_handle local) const {
+    assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton));
+    return degree_[local.vertex];
+  }
+
+  //@}
+
+  /////////////////////////////////////////////////////////////////////////////
+  /** @name Edges operations
+   */
+  //@{
+ public:
+  /**
+   * @brief return an edge handle if the two vertices forms
+   * an edge in the complex
+   */
+  boost::optional<Edge_handle> operator[](
+                                          const std::pair<Vertex_handle, Vertex_handle>& ab) const {
+    boost::optional<Edge_handle> res;
+    std::pair<Edge_handle, bool> edge_pair(
+                                           boost::edge(ab.first.vertex, ab.second.vertex, skeleton));
+    if (edge_pair.second)
+      res = edge_pair.first;
+    return res;
+  }
+
+  /**
+   * @brief returns the stored node associated to an edge
+   */
+  Graph_edge& operator[](Edge_handle edge_handle) {
+    return skeleton[edge_handle];
+  }
+
+  /**
+   * @brief returns the stored node associated to an edge
+   */
+  const Graph_edge& operator[](Edge_handle edge_handle) const {
+    return skeleton[edge_handle];
+  }
+
+  /**
+   * @brief returns the first vertex of an edge
+   * @details it assumes that the edge is present in the complex
+   */
+  Vertex_handle first_vertex(Edge_handle edge_handle) const {
+    return static_cast<Vertex_handle> (source(edge_handle, skeleton));
+  }
+
+  /**
+   * @brief returns the first vertex of an edge
+   * @details it assumes that the edge is present in the complex
+   */
+  Vertex_handle second_vertex(Edge_handle edge_handle) const {
+    return static_cast<Vertex_handle> (target(edge_handle, skeleton));
+  }
+
+  /**
+   * @brief returns the simplex made with the two vertices of the edge
+   * @details it assumes that the edge is present in the complex
+
+   */
+  Simplex_handle get_vertices(Edge_handle edge_handle) const {
+    auto edge((*this)[edge_handle]);
+    return Simplex_handle((*this)[edge.first()], (*this)[edge.second()]);
+  }
+
+  /**
+   * @brief Adds an edge between vertices a and b and all its cofaces.
+   */
+  Edge_handle add_edge(Vertex_handle a, Vertex_handle b) {
+    assert(contains_vertex(a) && contains_vertex(b));
+    assert(a != b);
+
+    auto edge_handle((*this)[std::make_pair(a, b)]);
+    // std::pair<Edge_handle,bool> pair_descr_bool = (*this)[std::make_pair(a,b)];
+    // Edge_handle edge_descr;
+    // bool edge_present = pair_descr_bool.second;
+    if (!edge_handle) {
+      edge_handle = boost::add_edge(a.vertex, b.vertex, skeleton).first;
+      (*this)[*edge_handle].setId(get_id(a), get_id(b));
+      degree_[a.vertex]++;
+      degree_[b.vertex]++;
+      if (visitor)
+        visitor->on_add_edge(a, b);
+    }
+    return *edge_handle;
+  }
+
+  /**
+   * @brief Adds all edges and their cofaces of a simplex to the simplicial complex.
+   */
+  void add_edges(const Simplex_handle & sigma) {
+    Simplex_handle_iterator i, j;
+    for (i = sigma.begin(); i != sigma.end(); ++i)
+      for (j = i, j++; j != sigma.end(); ++j)
+        add_edge(*i, *j);
+  }
+
+  /**
+   * @brief Removes an edge from the simplicial complex and all its cofaces.
+   * @details returns the former Edge_handle representing the edge
+   */
+  virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b) {
+    bool found;
+    Edge_handle edge;
+    tie(edge, found) = boost::edge(a.vertex, b.vertex, skeleton);
+    if (found) {
+      if (visitor)
+        visitor->on_remove_edge(a, b);
+      boost::remove_edge(a.vertex, b.vertex, skeleton);
+      degree_[a.vertex]--;
+      degree_[b.vertex]--;
+    }
+    return edge;
+  }
+
+  /**
+   * @brief Removes edge and its cofaces from the simplicial complex.
+   */
+  void remove_edge(Edge_handle edge) {
+    assert(contains_vertex(first_vertex(edge)));
+    assert(contains_vertex(second_vertex(edge)));
+    remove_edge(first_vertex(edge), second_vertex(edge));
+  }
+
+  /**
+   * @brief The complex is reduced to its set of vertices.
+   * All the edges and blockers are removed.
+   */
+  void keep_only_vertices() {
+    remove_blockers();
+
+    for (auto u : vertex_range()) {
+      while (this->degree(u) > 0) {
+        Vertex_handle v(*(adjacent_vertices(u.vertex, this->skeleton).first));
+        this->remove_edge(u, v);
+      }
+    }
+  }
+
+  /**
+   * @return true iff the simplicial complex contains an edge between
+   * vertices a and b
+   */
+  bool contains_edge(Vertex_handle a, Vertex_handle b) const {
+    // if (a.vertex<0 || b.vertex <0) return false;
+    return boost::edge(a.vertex, b.vertex, skeleton).second;
+  }
+
+  /**
+   * @return true iff the simplicial complex contains all vertices
+   * and all edges of simplex sigma
+   */
+  bool contains_edges(const Simplex_handle & sigma) const {
+    for (auto i = sigma.begin(); i != sigma.end(); ++i) {
+      if (!contains_vertex(*i))
+        return false;
+      for (auto j = i; ++j != sigma.end();) {
+        if (!contains_edge(*i, *j))
+          return false;
+      }
+    }
+    return true;
+  }
+  //@}
+
+  /////////////////////////////////////////////////////////////////////////////
+  /** @name Blockers operations
+   */
+  //@{
+
+  /**
+   * @brief Adds the simplex to the set of blockers and
+   * returns a Blocker_handle toward it if was not present before and 0 otherwise.
+   */
+  Blocker_handle add_blocker(const Simplex_handle& blocker) {
+    assert(blocker.dimension() > 1);
+    if (contains_blocker(blocker)) {
+      // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl;
+      return 0;
+    } else {
+      if (visitor)
+        visitor->on_add_blocker(blocker);
+      Blocker_handle blocker_pt = new Simplex_handle(blocker);
+      num_blockers_++;
+      auto vertex = blocker_pt->begin();
+      while (vertex != blocker_pt->end()) {
+        blocker_map_.insert(BlockerPair(*vertex, blocker_pt));
+        ++vertex;
+      }
+      return blocker_pt;
+    }
+  }
+
+ protected:
+  /**
+   * @brief Adds the simplex to the set of blockers
+   */
+  void add_blocker(Blocker_handle blocker) {
+    if (contains_blocker(*blocker)) {
+      // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl;
+      return;
+    } else {
+      if (visitor)
+        visitor->on_add_blocker(*blocker);
+      num_blockers_++;
+      auto vertex = blocker->begin();
+      while (vertex != blocker->end()) {
+        blocker_map_.insert(BlockerPair(*vertex, blocker));
+        ++vertex;
+      }
+    }
+  }
+
+ protected:
+  /**
+   * Removes sigma from the blocker map of vertex v
+   */
+  void remove_blocker(const Blocker_handle sigma, Vertex_handle v) {
+    Complex_blocker_around_vertex_iterator blocker;
+    for (blocker = blocker_range(v).begin(); blocker != blocker_range(v).end();
+         ++blocker) {
+      if (*blocker == sigma)
+        break;
+    }
+    if (*blocker != sigma) {
+      std::cerr
+          << "bug ((*blocker).second == sigma) ie try to remove a blocker not present\n";
+      assert(false);
+    } else {
+      blocker_map_.erase(blocker.current_position());
+    }
+  }
+
+ public:
+  /**
+   * @brief Removes the simplex from the set of blockers.
+   * @remark sigma has to belongs to the set of blockers
+   */
+  void remove_blocker(const Blocker_handle sigma) {
+    for (auto vertex : *sigma)
+      remove_blocker(sigma, vertex);
+    num_blockers_--;
+  }
+
+  /**
+   * @brief Remove all blockers, in other words, it expand the simplicial
+   * complex to the smallest flag complex that contains it.
+   */
+  void remove_blockers() {
+    // Desallocate the blockers
+    while (!blocker_map_.empty()) {
+      delete_blocker(blocker_map_.begin()->second);
+    }
+    num_blockers_ = 0;
+    blocker_map_.clear();
+  }
+
+ protected:
+  /**
+   * Removes the simplex sigma from the set of blockers.
+   * sigma has to belongs to the set of blockers
+   *
+   * @remark contrarily to delete_blockers does not call the destructor
+   */
+  void remove_blocker(const Simplex_handle& sigma) {
+    assert(contains_blocker(sigma));
+    for (auto vertex : sigma)
+      remove_blocker(sigma, vertex);
+    num_blockers_--;
+  }
+
+ public:
+  /**
+   * Removes the simplex s from the set of blockers
+   * and desallocate s.
+   */
+  void delete_blocker(Blocker_handle sigma) {
+    if (visitor)
+      visitor->on_delete_blocker(sigma);
+    remove_blocker(sigma);
+    delete sigma;
+  }
+
+  /**
+   * @return true iff s is a blocker of the simplicial complex
+   */
+  bool contains_blocker(const Blocker_handle s) const {
+    if (s->dimension() < 2)
+      return false;
+
+    Vertex_handle a = s->first_vertex();
+
+    for (const auto blocker : const_blocker_range(a)) {
+      if (s == *blocker)
+        return true;
+    }
+    return false;
+  }
+
+  /**
+   * @return true iff s is a blocker of the simplicial complex
+   */
+  bool contains_blocker(const Simplex_handle & s) const {
+    if (s.dimension() < 2)
+      return false;
+
+    Vertex_handle a = s.first_vertex();
+
+    for (auto blocker : const_blocker_range(a)) {
+      if (s == *blocker)
+        return true;
+    }
+    return false;
+  }
+
+ private:
+  /**
+   * @return true iff a blocker of the simplicial complex
+   * is a face of sigma.
+   */
+  bool blocks(const Simplex_handle & sigma) const {
+    for (auto s : sigma)
+      for (auto blocker : const_blocker_range(s))
+        if (sigma.contains(*blocker))
+          return true;
+    return false;
+  }
+
+  //@}
+
+ protected:
+  /**
+   * @details Adds to simplex the neighbours of v e.g. \f$ n \leftarrow n \cup N(v) \f$.
+   * If keep_only_superior is true then only vertices that are greater than v are added.
+   */
+  virtual void add_neighbours(Vertex_handle v, Simplex_handle & n,
+                              bool keep_only_superior = false) const {
+    boost_adjacency_iterator ai, ai_end;
+    for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end;
+         ++ai) {
+      if (keep_only_superior) {
+        if (*ai > v.vertex) {
+          n.add_vertex(Vertex_handle(*ai));
+        }
+      } else {
+        n.add_vertex(Vertex_handle(*ai));
+      }
+    }
+  }
+
+  /**
+   * @details Add to simplex res all vertices which are
+   * neighbours of alpha: ie \f$ res \leftarrow res \cup N(alpha) \f$.
+   *
+   * If 'keep_only_superior' is true then only vertices that are greater than alpha are added.
+   * todo revoir
+   *
+   */
+  virtual void add_neighbours(const Simplex_handle &alpha, Simplex_handle & res,
+                              bool keep_only_superior = false) const {
+    res.clear();
+    auto alpha_vertex = alpha.begin();
+    add_neighbours(*alpha_vertex, res, keep_only_superior);
+    for (alpha_vertex = (alpha.begin())++; alpha_vertex != alpha.end();
+         ++alpha_vertex)
+      keep_neighbours(*alpha_vertex, res, keep_only_superior);
+  }
+
+  /**
+   * @details Remove from simplex n all vertices which are
+   * not neighbours of v e.g. \f$ res \leftarrow res \cap N(v) \f$.
+   * If 'keep_only_superior' is true then only vertices that are greater than v are keeped.
+   */
+  virtual void keep_neighbours(Vertex_handle v, Simplex_handle& res,
+                               bool keep_only_superior = false) const {
+    Simplex_handle nv;
+    add_neighbours(v, nv, keep_only_superior);
+    res.intersection(nv);
+  }
+
+  /**
+   * @details Remove from simplex all vertices which are
+   * neighbours of v eg \f$ res \leftarrow res \setminus N(v) \f$.
+   * If 'keep_only_superior' is true then only vertices that are greater than v are added.
+   */
+  virtual void remove_neighbours(Vertex_handle v, Simplex_handle & res,
+                                 bool keep_only_superior = false) const {
+    Simplex_handle nv;
+    add_neighbours(v, nv, keep_only_superior);
+    res.difference(nv);
+  }
+
+ public:
+  typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex> Link_complex;
+
+  /**
+   * Constructs the link of 'simplex' with points coordinates.
+   */
+  Link_complex link(Vertex_handle v) const {
+    return Link_complex(*this, Simplex_handle(v));
+  }
+
+  /**
+   * Constructs the link of 'simplex' with points coordinates.
+   */
+  Link_complex link(Edge_handle edge) const {
+    return Link_complex(*this, edge);
+  }
+
+  /**
+   * Constructs the link of 'simplex' with points coordinates.
+   */
+  Link_complex link(const Simplex_handle& simplex) const {
+    return Link_complex(*this, simplex);
+  }
+
+  /**
+   * @brief Compute the local vertices of 's' in the current complex
+   * If one of them is not present in the complex then the return value is uninitialized.
+   *
+   *
+   */
+  // xxx rename get_address et place un using dans sub_complex
+
+  boost::optional<Simplex_handle> get_simplex_address(
+                                                      const Root_simplex_handle& s) const {
+    boost::optional<Simplex_handle> res;
+
+    Simplex_handle s_address;
+    // Root_simplex_const_iterator i;
+    for (auto i = s.begin(); i != s.end(); ++i) {
+      boost::optional<Vertex_handle> address = get_address(*i);
+      if (!address)
+        return res;
+      else
+        s_address.add_vertex(*address);
+    }
+    res = s_address;
+    return res;
+  }
+
+  /**
+   * @brief returns a simplex with vertices which are the id of vertices of the
+   * argument.
+   */
+  Root_simplex_handle get_id(const Simplex_handle& local_simplex) const {
+    Root_simplex_handle global_simplex;
+    for (auto x = local_simplex.begin(); x != local_simplex.end(); ++x) {
+      global_simplex.add_vertex(get_id(*x));
+    }
+    return global_simplex;
+  }
+
+  /**
+   * @brief returns true iff the simplex s belongs to the simplicial
+   * complex.
+   */
+  virtual bool contains(const Simplex_handle & s) const {
+    if (s.dimension() == -1) {
+      return false;
+    } else if (s.dimension() == 0) {
+      return contains_vertex(s.first_vertex());
+    } else {
+      return (contains_edges(s) && !blocks(s));
+    }
+  }
+
+  /*
+   * @brief returnrs true iff the complex is empty.
+   */
+  bool empty() const {
+    return num_vertices() == 0;
+  }
+
+  /*
+   * @brief returns the number of vertices in the complex.
+   */
+  int num_vertices() const {
+    // remark boost::num_vertices(skeleton) counts deactivated vertices
+    return num_vertices_;
+  }
+
+  /*
+   * @brief returns the number of edges in the complex.
+   * @details currently in O(n)
+   */
+  // todo cache the value
+
+  int num_edges() const {
+    return boost::num_edges(skeleton);
+  }
+
+  int num_triangles() const {
+    auto triangles = triangle_range();
+    return std::distance(triangles.begin(), triangles.end());
+  }
+
+  /*
+   * @brief returns the number of blockers in the complex.
+   */
+  int num_blockers() const {
+    return num_blockers_;
+  }
+
+  /*
+   * @brief returns true iff the graph of the 1-skeleton of the complex is complete.
+   */
+  bool complete() const {
+    return (num_vertices() * (num_vertices() - 1)) / 2 == num_edges();
+  }
+
+  /**
+   * @brief returns the number of connected components in the graph of the 1-skeleton.
+   */
+  int num_connected_components() const {
+    int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices();
+    std::vector<int> component(skeleton.vertex_set().size());
+    return boost::connected_components(this->skeleton, &component[0])
+        - num_vert_collapsed;
+  }
+
+  /**
+   * @brief %Test if the complex is a cone.
+   * @details Runs in O(n) where n is the number of vertices.
+   */
+  bool is_cone() const {
+    if (num_vertices() == 0)
+      return false;
+    if (num_vertices() == 1)
+      return true;
+    for (auto vi : vertex_range()) {
+      // xxx todo faire une methode bool is_in_blocker(Vertex_handle)
+      if (blocker_map_.find(vi) == blocker_map_.end()) {
+        // no blocker passes through the vertex, we just need to
+        // check if the current vertex is linked to all others vertices of the complex
+        if (degree_[vi.vertex] == num_vertices() - 1)
+          return true;
+      }
+    }
+    return false;
+  }
+
+  //@}
+  /** @Simplification operations
+   */
+  //@{
+
+  /**
+   * Returns true iff the blocker 'sigma' is popable.
+   * To define popable, let us call 'L' the complex that
+   * consists in the current complex without the blocker 'sigma'.
+   * A blocker 'sigma' is then "popable" if the link of 'sigma'
+   * in L is reducible.
+   *
+   */
+  bool is_popable_blocker(Blocker_handle sigma) const;
+
+  /**
+   * Removes all the popable blockers of the complex and delete them.
+   * @returns the number of popable blockers deleted
+   */
+  void remove_popable_blockers();
+
+  /**
+   * Removes all the popable blockers of the complex passing through v and delete them.
+   */
+  void remove_popable_blockers(Vertex_handle v);
+
+  /**
+   * @brief Removes all the popable blockers of the complex passing through v and delete them.
+   * Also remove popable blockers in the neighborhood if they became popable.
+   *
+   */
+  void remove_all_popable_blockers(Vertex_handle v);
+
+  /**
+   * Remove the star of the vertex 'v'
+   */
+  void remove_star(Vertex_handle v);
+
+ private:
+  /**
+   * after removing the star of a simplex, blockers sigma that contains this simplex must be removed.
+   * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added
+   * whenever the dimension of tau is at least 2.
+   */
+  void update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed);
+
+ public:
+  /**
+   * Remove the star of the edge connecting vertices a and b.
+   * @returns the number of blocker that have been removed
+   */
+  void remove_star(Vertex_handle a, Vertex_handle b);
+
+  /**
+   * Remove the star of the edge 'e'.
+   */
+  void remove_star(Edge_handle e);
+
+  /**
+   * Remove the star of the simplex 'sigma' which needs to belong to the complex
+   */
+  void remove_star(const Simplex_handle& sigma);
+
+  /**
+   * @brief add a maximal simplex plus all its cofaces.
+   * @details the simplex must have dimension greater than one (otherwise use add_vertex or add_edge).
+   */
+  void add_simplex(const Simplex_handle& sigma);
+
+ private:
+  /**
+   * remove all blockers that contains sigma
+   */
+  void remove_blocker_containing_simplex(const Simplex_handle& sigma);
+
+  /**
+   * remove all blockers that contains sigma
+   */
+  void remove_blocker_include_in_simplex(const Simplex_handle& sigma);
+
+ public:
+  enum simplifiable_status {
+    NOT_HOMOTOPY_EQ, MAYBE_HOMOTOPY_EQ, HOMOTOPY_EQ
+  };
+
+  simplifiable_status is_remove_star_homotopy_preserving(const Simplex_handle& simplex) {
+    // todo write a virtual method 'link' in Skeleton_blocker_complex which will be overloaded by the current one of
+    // Skeleton_blocker_geometric_complex
+    // then call it there to build the link and return the value of link.is_contractible()
+    return MAYBE_HOMOTOPY_EQ;
+  }
+
+  enum contractible_status {
+    NOT_CONTRACTIBLE, MAYBE_CONTRACTIBLE, CONTRACTIBLE
+  };
+
+  /**
+   * @brief %Test if the complex is reducible using a strategy defined in the class
+   * (by default it tests if the complex is a cone)
+   * @details Note that NO could be returned if some invariant ensures that the complex
+   * is not a point (for instance if the euler characteristic is different from 1).
+   * This function will surely have to return MAYBE in some case because the
+   * associated problem is undecidable but it in practice, it can often
+   * be solved with the help of geometry.
+   */
+  virtual contractible_status is_contractible() const {
+    if (this->is_cone()) {
+      return CONTRACTIBLE;
+    } else {
+      return MAYBE_CONTRACTIBLE;
+    }
+  }
+  //@}
+
+  /** @Edge contraction operations
+   */
+  //@{
+
+  /**
+   * @return If ignore_popable_blockers is true
+   * then the result is true iff the link condition at edge ab is satisfied
+   * or equivalently iff no blocker contains ab.
+   * If ignore_popable_blockers is false then the
+   * result is true iff all blocker containing ab are popable.
+   */
+  bool link_condition(Vertex_handle a, Vertex_handle b, bool ignore_popable_blockers = false) const {
+    for (auto blocker : this->const_blocker_range(a))
+      if (blocker->contains(b)) {
+        // false if ignore_popable_blockers is false
+        // otherwise the blocker has to be popable
+        return ignore_popable_blockers && is_popable_blocker(blocker);
+      }
+    return true;
+  }
+
+  /**
+   * @return If ignore_popable_blockers is true
+   * then the result is true iff the link condition at edge ab is satisfied
+   * or equivalently iff no blocker contains ab.
+   * If ignore_popable_blockers is false then the
+   * result is true iff all blocker containing ab are popable.
+   */
+  bool link_condition(Edge_handle e, bool ignore_popable_blockers = false) const {
+    const Graph_edge& edge = (*this)[e];
+    assert(this->get_address(edge.first()));
+    assert(this->get_address(edge.second()));
+    Vertex_handle a(*this->get_address(edge.first()));
+    Vertex_handle b(*this->get_address(edge.second()));
+    return link_condition(a, b, ignore_popable_blockers);
+  }
+
+ protected:
+  /**
+   * Compute simplices beta such that a.beta is an order 0 blocker
+   * that may be used to construct a new blocker after contracting ab.
+   * It requires that the link condition is satisfied.
+   */
+  void tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const;
+
+ private:
+  /**
+   * @brief "Replace" the edge 'bx' by the edge 'ax'.
+   * Assume that the edge 'bx' was present whereas 'ax' was not.
+   * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'.
+   * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added
+   * and just before edge 'bx' had been removed. That way, it can
+   * eventually access to information of 'ax'.
+   */
+  void swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x);
+
+ private:
+  /**
+   * @brief removes all blockers passing through the edge 'ab'
+   */
+  void delete_blockers_around_vertex(Vertex_handle v);
+
+  /**
+   * @brief removes all blockers passing through the edge 'ab'
+   */
+  void delete_blockers_around_edge(Vertex_handle a, Vertex_handle b);
+
+ public:
+  /**
+   * Contracts the edge.
+   * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
+   * it removes first all blockers passing through 'ab'
+   */
+  void contract_edge(Edge_handle edge) {
+    contract_edge(this->first_vertex(edge), this->second_vertex(edge));
+  }
+
+  /**
+   * Contracts the edge connecting vertices a and b.
+   * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
+   * it removes first all blockers passing through 'ab'
+   */
+  void contract_edge(Vertex_handle a, Vertex_handle b);
+
+ private:
+  void get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b,
+                                                  std::set<Simplex_handle>& blockers_to_add);
+  /**
+   * delete all blockers that passes through a or b
+   */
+  void delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b);
+  void update_edges_after_contraction(Vertex_handle a, Vertex_handle b);
+  void notify_changed_edges(Vertex_handle a);
+  //@}
+
+ public:
+  /////////////////////////////////////////////////////////////////////////////
+  /** @name Vertex iterators
+   */
+  //@{
+  typedef Vertex_iterator<Skeleton_blocker_complex> Complex_vertex_iterator;
+
+  //
+  // Range over the vertices of the simplicial complex.
+  // Methods .begin() and .end() return a Complex_vertex_iterator.
+  //
+  typedef boost::iterator_range<Complex_vertex_iterator> Complex_vertex_range;
+
+  /**
+   * @brief Returns a Complex_vertex_range over all vertices of the complex
+   */
+  Complex_vertex_range vertex_range() const {
+    auto begin = Complex_vertex_iterator(this);
+    auto end = Complex_vertex_iterator(this, 0);
+    return Complex_vertex_range(begin, end);
+  }
+
+  typedef Neighbors_vertices_iterator<Skeleton_blocker_complex> Complex_neighbors_vertices_iterator;
+
+
+  typedef boost::iterator_range<Complex_neighbors_vertices_iterator> Complex_neighbors_vertices_range;
+
+  /**
+   * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v
+   */
+  Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const {
+    auto begin = Complex_neighbors_vertices_iterator(this, v);
+    auto end = Complex_neighbors_vertices_iterator(this, v, 0);
+    return Complex_neighbors_vertices_range(begin, end);
+  }
+
+  //@}
+
+  /** @name Edge iterators
+   */
+  //@{
+
+  typedef Edge_iterator<Skeleton_blocker_complex> Complex_edge_iterator;
+
+
+  typedef boost::iterator_range<Complex_edge_iterator> Complex_edge_range;
+
+  /**
+   * @brief Returns a Complex_edge_range over all edges of the simplicial complex
+   */
+  Complex_edge_range edge_range() const {
+    auto begin = Complex_edge_iterator(this);
+    auto end = Complex_edge_iterator(this, 0);
+    return Complex_edge_range(begin, end);
+  }
+
+
+  typedef Edge_around_vertex_iterator<Skeleton_blocker_complex> Complex_edge_around_vertex_iterator;
+
+
+  typedef boost::iterator_range <Complex_edge_around_vertex_iterator> Complex_edge_around_vertex_range;
+
+  /**
+   * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes
+   * through 'v'
+   */
+  Complex_edge_around_vertex_range edge_range(Vertex_handle v) const {
+    auto begin = Complex_edge_around_vertex_iterator(this, v);
+    auto end = Complex_edge_around_vertex_iterator(this, v, 0);
+    return Complex_edge_around_vertex_range(begin, end);
+  }
+
+  //@}
+
+  /** @name Triangles iterators
+   */
+  //@{
+ private:
+  typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link;
+  typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link;
+
+ public:
+  typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex, Superior_link>
+  Superior_triangle_around_vertex_iterator;
+  typedef boost::iterator_range < Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link> >
+  Complex_triangle_around_vertex_range;
+
+  /**
+   * @brief Range over triangles around a vertex of the simplicial complex.
+   * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
+   *
+   */
+  Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const {
+    auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v);
+    auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v, 0);
+    return Complex_triangle_around_vertex_range(begin, end);
+  }
+
+  typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range;
+  typedef Triangle_iterator<Skeleton_blocker_complex> Complex_triangle_iterator;
+
+  /**
+   * @brief Range over triangles of the simplicial complex.
+   * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
+   *
+   */
+  Complex_triangle_range triangle_range() const {
+    auto end = Triangle_iterator<Skeleton_blocker_complex>(this, 0);
+    if (empty()) {
+      return Complex_triangle_range(end, end);
+    } else {
+      auto begin = Triangle_iterator<Skeleton_blocker_complex>(this);
+      return Complex_triangle_range(begin, end);
+    }
+  }
+
+  //@}
+
+  /** @name Simplices iterators
+   */
+  //@{
+  typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex, Link> Complex_simplex_around_vertex_iterator;
+
+  /**
+   * @brief Range over the simplices of the simplicial complex around a vertex.
+   * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator.
+   */
+  typedef boost::iterator_range < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range;
+
+  /**
+   * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex
+   */
+  Complex_simplex_around_vertex_range simplex_range(Vertex_handle v) const {
+    assert(contains_vertex(v));
+    return Complex_simplex_around_vertex_range(
+                                               Complex_simplex_around_vertex_iterator(this, v),
+                                               Complex_simplex_around_vertex_iterator(this, v, true));
+  }
+
+  // typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator;
+  typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator;
+
+  typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range;
+
+  /**
+   * @brief Returns a Complex_simplex_range over all the simplices of the complex
+   */
+  Complex_simplex_range simplex_range() const {
+    Complex_simplex_iterator end(this, true);
+    if (empty()) {
+      return Complex_simplex_range(end, end);
+    } else {
+      Complex_simplex_iterator begin(this);
+      return Complex_simplex_range(begin, end);
+    }
+  }
+
+  //@}
+
+  /** @name Blockers iterators
+   */
+  //@{
+ private:
+  /**
+   * @brief Iterator over the blockers adjacent to a vertex
+   */
+  typedef Blocker_iterator_around_vertex_internal<
+  typename std::multimap<Vertex_handle, Simplex_handle *>::iterator,
+  Blocker_handle>
+  Complex_blocker_around_vertex_iterator;
+
+  /**
+   * @brief Iterator over (constant) blockers adjacent to a vertex
+   */
+  typedef Blocker_iterator_around_vertex_internal<
+  typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator,
+  const Blocker_handle>
+  Const_complex_blocker_around_vertex_iterator;
+
+  typedef boost::iterator_range <Complex_blocker_around_vertex_iterator> Complex_blocker_around_vertex_range;
+  typedef boost::iterator_range <Const_complex_blocker_around_vertex_iterator> Const_complex_blocker_around_vertex_range;
+
+ public:
+  /**
+   * @brief Returns a range of the blockers of the complex passing through a vertex
+   */
+  Complex_blocker_around_vertex_range blocker_range(Vertex_handle v) {
+    auto begin = Complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v));
+    auto end = Complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v));
+    return Complex_blocker_around_vertex_range(begin, end);
+  }
+
+  /**
+   * @brief Returns a range of the blockers of the complex passing through a vertex
+   */
+  Const_complex_blocker_around_vertex_range const_blocker_range(Vertex_handle v) const {
+    auto begin = Const_complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v));
+    auto end = Const_complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v));
+    return Const_complex_blocker_around_vertex_range(begin, end);
+  }
+
+ private:
+  /**
+   * @brief Iterator over the blockers.
+   */
+  typedef Blocker_iterator_internal<
+  typename std::multimap<Vertex_handle, Simplex_handle *>::iterator,
+  Blocker_handle>
+  Complex_blocker_iterator;
+
+  /**
+   * @brief Iterator over the (constant) blockers.
+   */
+  typedef Blocker_iterator_internal<
+  typename std::multimap<Vertex_handle, Simplex_handle *>::const_iterator,
+  const Blocker_handle>
+  Const_complex_blocker_iterator;
+
+  typedef boost::iterator_range <Complex_blocker_iterator> Complex_blocker_range;
+  typedef boost::iterator_range <Const_complex_blocker_iterator> Const_complex_blocker_range;
+
+ public:
+  /**
+   * @brief Returns a range of the blockers of the complex
+   */
+  Complex_blocker_range blocker_range() {
+    auto begin = Complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end());
+    auto end = Complex_blocker_iterator(blocker_map_.end(), blocker_map_.end());
+    return Complex_blocker_range(begin, end);
+  }
+
+  /**
+   * @brief Returns a range of the blockers of the complex
+   */
+  Const_complex_blocker_range const_blocker_range() const {
+    auto begin = Const_complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end());
+    auto end = Const_complex_blocker_iterator(blocker_map_.end(), blocker_map_.end());
+    return Const_complex_blocker_range(begin, end);
+  }
+
+  //@}
+
+  /////////////////////////////////////////////////////////////////////////////
+  /** @name Print and IO methods
+   */
+  //@{
+ public:
+  std::string to_string() const {
+    std::ostringstream stream;
+    stream << num_vertices() << " vertices:\n" << vertices_to_string() << std::endl;
+    stream << num_edges() << " edges:\n" << edges_to_string() << std::endl;
+    stream << num_blockers() << " blockers:\n" << blockers_to_string() << std::endl;
+    return stream.str();
+  }
+
+  std::string vertices_to_string() const {
+    std::ostringstream stream;
+    for (auto vertex : vertex_range()) {
+      stream << "{" << (*this)[vertex].get_id() << "} ";
+    }
+    stream << std::endl;
+    return stream.str();
+  }
+
+  std::string edges_to_string() const {
+    std::ostringstream stream;
+    for (auto edge : edge_range())
+      stream << "{" << (*this)[edge].first() << "," << (*this)[edge].second() << "} ";
+    stream << std::endl;
+    return stream.str();
+  }
+
+  std::string blockers_to_string() const {
+    std::ostringstream stream;
+
+    for (auto b : const_blocker_range())
+      stream << *b << std::endl;
+    return stream.str();
+  }
+  //@}
 };
 
-
-
 /**
  * build a simplicial complex from a collection
  * of top faces.
  * return the total number of simplices
  */
 template<typename Complex, typename SimplexHandleIterator>
-Complex make_complex_from_top_faces(SimplexHandleIterator simplex_begin, SimplexHandleIterator simplex_end, bool is_flag_complex = false) {
-	typedef typename Complex::Simplex_handle Simplex_handle;
-	std::vector<Simplex_handle> simplices;
-	for (auto top_face = simplex_begin; top_face != simplex_end; ++top_face) {
-		auto subfaces_topface = subfaces(*top_face);
-		simplices.insert(simplices.end(),subfaces_topface.begin(),subfaces_topface.end());
-	}
-	return Complex(simplices.begin(),simplices.end(),is_flag_complex);
+Complex make_complex_from_top_faces(SimplexHandleIterator simplex_begin, SimplexHandleIterator simplex_end, 
+                                    bool is_flag_complex = false) {
+  typedef typename Complex::Simplex_handle Simplex_handle;
+  std::vector<Simplex_handle> simplices;
+  for (auto top_face = simplex_begin; top_face != simplex_end; ++top_face) {
+    auto subfaces_topface = subfaces(*top_face);
+    simplices.insert(simplices.end(), subfaces_topface.begin(), subfaces_topface.end());
+  }
+  return Complex(simplices.begin(), simplices.end(), is_flag_complex);
 }
 
 }  // namespace skbl
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
index 69de1085..437482cb 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h
@@ -37,195 +37,186 @@ namespace skbl {
  */
 template<typename SkeletonBlockerGeometricDS>
 class Skeleton_blocker_geometric_complex :
-		public Skeleton_blocker_complex<SkeletonBlockerGeometricDS> {
-		public:
-	typedef typename SkeletonBlockerGeometricDS::GT GT;
-
-	typedef Skeleton_blocker_complex<SkeletonBlockerGeometricDS> SimplifiableSkeletonblocker;
-
-	typedef typename SimplifiableSkeletonblocker::Vertex_handle Vertex_handle;
-	typedef typename SimplifiableSkeletonblocker::Root_vertex_handle Root_vertex_handle;
-	typedef typename SimplifiableSkeletonblocker::Edge_handle Edge_handle;
-	typedef typename SimplifiableSkeletonblocker::Simplex_handle Simplex_handle;
-
-	typedef typename SimplifiableSkeletonblocker::Graph_vertex Graph_vertex;
-
-	typedef typename SkeletonBlockerGeometricDS::Point Point;
-
-
-
-	Skeleton_blocker_geometric_complex(){
-	}
-
-
-
-	/**
-	 * constructor given a list of points
-	 */
-	template<typename PointIterator>
-	explicit Skeleton_blocker_geometric_complex(int num_vertices,PointIterator begin,PointIterator end){
-		for(auto point = begin; point != end; ++point)
-			add_vertex(*point);
-	}
-
-
-	/**
-	 * @brief Constructor with a list of simplices.
-	 * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not).
-	 */
-	template<typename SimpleHandleOutputIterator,typename PointIterator>
-	Skeleton_blocker_geometric_complex(
-			SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end,
-			PointIterator points_begin,PointIterator points_end,
-			bool is_flag_complex = false):Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin,simplex_end,is_flag_complex){
-		unsigned current = 0;
-		for(auto point = points_begin; point != points_end; ++point)
-			(*this)[Vertex_handle(current++)].point() = Point(point->begin(),point->end());
-	}
-
-	template<typename SimpleHandleOutputIterator,typename PointIterator>
-	friend Skeleton_blocker_geometric_complex make_complex_from_top_faces(
-			SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end,
-			PointIterator points_begin,PointIterator points_end,
-			bool is_flag_complex = false){
-		Skeleton_blocker_geometric_complex complex;
-		unsigned current = 0;
-		complex=make_complex_from_top_faces<Skeleton_blocker_geometric_complex>(simplex_begin,simplex_end,is_flag_complex);
-		for(auto point = points_begin; point != points_end; ++point)
-			//			complex.point(Vertex_handle(current++)) = Point(point->begin(),point->end());
-			complex.point(Vertex_handle(current++)) = Point(*point);
-		return complex;
-	}
-
-
-	/**
-	 * @brief Constructor with a list of simplices.
-	 * Points of every vertex are the point constructed with default constructor.
-	 * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not).
-	 */
-	template<typename SimpleHandleOutputIterator>
-	Skeleton_blocker_geometric_complex(
-			SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end,
-			bool is_flag_complex = false):Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin,simplex_end,is_flag_complex){
-	}
-
-
-	/**
-	 * @brief Add a vertex to the complex with a default constructed associated point.
-	 */
-	Vertex_handle add_vertex() {
-		return SimplifiableSkeletonblocker::add_vertex();
-	}
-
-	/**
-	 * @brief Add a vertex to the complex with its associated point.
-	 */
-	Vertex_handle add_vertex(const Point& point) {
-		Vertex_handle ad = SimplifiableSkeletonblocker::add_vertex();
-		(*this)[ad].point() = point;
-		return ad;
-	}
-
-	/**
-	 * @brief Returns the Point associated to the vertex v.
-	 */
-	const Point& point(Vertex_handle v) const {
-		assert(this->contains_vertex(v));
-		return (*this)[v].point();
-	}
-
-	/**
-	 * @brief Returns the Point associated to the vertex v.
-	 */
-	Point& point(Vertex_handle v) {
-		assert(this->contains_vertex(v));
-		return (*this)[v].point();
-	}
-
-	const Point& point(Root_vertex_handle global_v) const {
-		Vertex_handle local_v((*this)[global_v]);
-		assert(this->contains_vertex(local_v));
-		return (*this)[local_v].point();
-	}
-
-	Point& point(Root_vertex_handle global_v) {
-		Vertex_handle local_v((*this)[global_v]);
-		assert(this->contains_vertex(local_v));
-		return (*this)[local_v].point();
-	}
-
-	typedef Skeleton_blocker_link_complex<Skeleton_blocker_geometric_complex> Geometric_link;
-
-	/**
-	 * Constructs the link of 'simplex' with points coordinates.
-	 */
-	Geometric_link link(Vertex_handle v) const {
-		Geometric_link link(*this, Simplex_handle(v));
-		// we now add the point info
-		add_points_to_link(link);
-		return link;
-	}
-
-	/**
-	 * Constructs the link of 'simplex' with points coordinates.
-	 */
-	Geometric_link link(const Simplex_handle& simplex) const {
-		Geometric_link link(*this, simplex);
-		// we now add the point info
-		add_points_to_link(link);
-		return link;
-	}
-
-	/**
-	 * Constructs the link of 'simplex' with points coordinates.
-	 */
-	Geometric_link link(Edge_handle edge) const {
-		Geometric_link link(*this, edge);
-		// we now add the point info
-		add_points_to_link(link);
-		return link;
-	}
-
-	typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerGeometricDS>> Abstract_link;
-
-	/**
-	 * Constructs the abstract link of v (without points coordinates).
-	 */
-	Abstract_link abstract_link(Vertex_handle v) const {
-		return Abstract_link(*this, Simplex_handle(v));
-	}
-
-	/**
-	 * Constructs the link of 'simplex' with points coordinates.
-	 */
-	Abstract_link abstract_link(const Simplex_handle& simplex) const {
-		return Abstract_link(*this, simplex);
-	}
-
-	/**
-	 * Constructs the link of 'simplex' with points coordinates.
-	 */
-	Abstract_link abstract_link(Edge_handle edge) const {
-		return Abstract_link(*this, edge);
-	}
-
-
-
-		private:
-
-	void add_points_to_link(Geometric_link& link) const {
-		for (Vertex_handle v : link.vertex_range()) {
-			Root_vertex_handle v_root(link.get_id(v));
-			link.point(v) = (*this).point(v_root);
-		}
-	}
-
-
-
+public Skeleton_blocker_complex<SkeletonBlockerGeometricDS> {
+ public:
+  typedef typename SkeletonBlockerGeometricDS::GT GT;
+
+  typedef Skeleton_blocker_complex<SkeletonBlockerGeometricDS> SimplifiableSkeletonblocker;
+
+  typedef typename SimplifiableSkeletonblocker::Vertex_handle Vertex_handle;
+  typedef typename SimplifiableSkeletonblocker::Root_vertex_handle Root_vertex_handle;
+  typedef typename SimplifiableSkeletonblocker::Edge_handle Edge_handle;
+  typedef typename SimplifiableSkeletonblocker::Simplex_handle Simplex_handle;
+
+  typedef typename SimplifiableSkeletonblocker::Graph_vertex Graph_vertex;
+
+  typedef typename SkeletonBlockerGeometricDS::Point Point;
+
+  Skeleton_blocker_geometric_complex() { }
+
+  /**
+   * constructor given a list of points
+   */
+  template<typename PointIterator>
+  explicit Skeleton_blocker_geometric_complex(int num_vertices, PointIterator begin, PointIterator end) {
+    for (auto point = begin; point != end; ++point)
+      add_vertex(*point);
+  }
+
+  /**
+   * @brief Constructor with a list of simplices.
+   * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be
+   * computed or not).
+   */
+  template<typename SimpleHandleOutputIterator, typename PointIterator>
+  Skeleton_blocker_geometric_complex(
+                                     SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end,
+                                     PointIterator points_begin, PointIterator points_end,
+                                     bool is_flag_complex = false)
+      : Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin, simplex_end, is_flag_complex) {
+    unsigned current = 0;
+    for (auto point = points_begin; point != points_end; ++point)
+      (*this)[Vertex_handle(current++)].point() = Point(point->begin(), point->end());
+  }
+
+  template<typename SimpleHandleOutputIterator, typename PointIterator>
+  friend Skeleton_blocker_geometric_complex make_complex_from_top_faces(
+                                                                        SimpleHandleOutputIterator simplex_begin,
+                                                                        SimpleHandleOutputIterator simplex_end,
+                                                                        PointIterator points_begin,
+                                                                        PointIterator points_end,
+                                                                        bool is_flag_complex = false) {
+    Skeleton_blocker_geometric_complex complex;
+    unsigned current = 0;
+    complex =
+        make_complex_from_top_faces<Skeleton_blocker_geometric_complex>(simplex_begin, simplex_end, is_flag_complex);
+    for (auto point = points_begin; point != points_end; ++point)
+      // complex.point(Vertex_handle(current++)) = Point(point->begin(),point->end());
+      complex.point(Vertex_handle(current++)) = Point(*point);
+    return complex;
+  }
+
+  /**
+   * @brief Constructor with a list of simplices.
+   * Points of every vertex are the point constructed with default constructor.
+   * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not).
+   */
+  template<typename SimpleHandleOutputIterator>
+  Skeleton_blocker_geometric_complex(
+                                     SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end,
+                                     bool is_flag_complex = false)
+      : Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin, simplex_end, is_flag_complex) { }
+
+  /**
+   * @brief Add a vertex to the complex with a default constructed associated point.
+   */
+  Vertex_handle add_vertex() {
+    return SimplifiableSkeletonblocker::add_vertex();
+  }
+
+  /**
+   * @brief Add a vertex to the complex with its associated point.
+   */
+  Vertex_handle add_vertex(const Point& point) {
+    Vertex_handle ad = SimplifiableSkeletonblocker::add_vertex();
+    (*this)[ad].point() = point;
+    return ad;
+  }
+
+  /**
+   * @brief Returns the Point associated to the vertex v.
+   */
+  const Point& point(Vertex_handle v) const {
+    assert(this->contains_vertex(v));
+    return (*this)[v].point();
+  }
+
+  /**
+   * @brief Returns the Point associated to the vertex v.
+   */
+  Point& point(Vertex_handle v) {
+    assert(this->contains_vertex(v));
+    return (*this)[v].point();
+  }
+
+  const Point& point(Root_vertex_handle global_v) const {
+    Vertex_handle local_v((*this)[global_v]);
+    assert(this->contains_vertex(local_v));
+    return (*this)[local_v].point();
+  }
+
+  Point& point(Root_vertex_handle global_v) {
+    Vertex_handle local_v((*this)[global_v]);
+    assert(this->contains_vertex(local_v));
+    return (*this)[local_v].point();
+  }
+
+  typedef Skeleton_blocker_link_complex<Skeleton_blocker_geometric_complex> Geometric_link;
+
+  /**
+   * Constructs the link of 'simplex' with points coordinates.
+   */
+  Geometric_link link(Vertex_handle v) const {
+    Geometric_link link(*this, Simplex_handle(v));
+    // we now add the point info
+    add_points_to_link(link);
+    return link;
+  }
+
+  /**
+   * Constructs the link of 'simplex' with points coordinates.
+   */
+  Geometric_link link(const Simplex_handle& simplex) const {
+    Geometric_link link(*this, simplex);
+    // we now add the point info
+    add_points_to_link(link);
+    return link;
+  }
+
+  /**
+   * Constructs the link of 'simplex' with points coordinates.
+   */
+  Geometric_link link(Edge_handle edge) const {
+    Geometric_link link(*this, edge);
+    // we now add the point info
+    add_points_to_link(link);
+    return link;
+  }
+
+  typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerGeometricDS>> Abstract_link;
+
+  /**
+   * Constructs the abstract link of v (without points coordinates).
+   */
+  Abstract_link abstract_link(Vertex_handle v) const {
+    return Abstract_link(*this, Simplex_handle(v));
+  }
+
+  /**
+   * Constructs the link of 'simplex' with points coordinates.
+   */
+  Abstract_link abstract_link(const Simplex_handle& simplex) const {
+    return Abstract_link(*this, simplex);
+  }
+
+  /**
+   * Constructs the link of 'simplex' with points coordinates.
+   */
+  Abstract_link abstract_link(Edge_handle edge) const {
+    return Abstract_link(*this, edge);
+  }
+
+ private:
+  void add_points_to_link(Geometric_link& link) const {
+    for (Vertex_handle v : link.vertex_range()) {
+      Root_vertex_handle v_root(link.get_id(v));
+      link.point(v) = (*this).point(v_root);
+    }
+  }
 };
 
-} // namespace skbl
+}  // namespace skbl
 
-} // namespace Gudhi
+}  // namespace Gudhi
 
 #endif  // SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
index 9eab7f1e..86a12d90 100644
--- a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
+++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h
@@ -32,316 +32,52 @@ namespace Gudhi {
 
 namespace skbl {
 
-///**
-// *  \brief Class that allows simplification operation on a simplicial complex represented
-// *  by a skeleton/blockers pair.
-// * \ingroup skbl
-// * @extends Skeleton_blocker_complex
-// */
-//template<typename SkeletonBlockerDS>
-//class Skeleton_blocker_complex : public Skeleton_blocker_complex<SkeletonBlockerDS> {
-//	template<class ComplexType> friend class Skeleton_blocker_sub_complex;
-//
-//public:
-//	typedef Skeleton_blocker_complex<SkeletonBlockerDS> SkeletonBlockerComplex;
-//
-//	typedef typename SkeletonBlockerComplex::Graph_edge Graph_edge;
-//
-//	typedef typename SkeletonBlockerComplex::boost_adjacency_iterator boost_adjacency_iterator;
-//	typedef typename SkeletonBlockerComplex::Edge_handle Edge_handle;
-//	typedef typename SkeletonBlockerComplex::boost_vertex_handle boost_vertex_handle;
-//	typedef typename SkeletonBlockerComplex::Vertex_handle Vertex_handle;
-//	typedef typename SkeletonBlockerComplex::Root_vertex_handle Root_vertex_handle;
-//	typedef typename SkeletonBlockerComplex::Simplex_handle Simplex_handle;
-//	typedef typename SkeletonBlockerComplex::Root_simplex_handle Root_simplex_handle;
-//	typedef typename SkeletonBlockerComplex::Blocker_handle Blocker_handle;
-//
-//
-//	typedef typename SkeletonBlockerComplex::Root_simplex_iterator Root_simplex_iterator;
-//	typedef typename SkeletonBlockerComplex::Simplex_handle_iterator Simplex_handle_iterator;
-//	typedef typename SkeletonBlockerComplex::BlockerMap BlockerMap;
-//	typedef typename SkeletonBlockerComplex::BlockerPair BlockerPair;
-//	typedef typename SkeletonBlockerComplex::BlockerMapIterator BlockerMapIterator;
-//	typedef typename SkeletonBlockerComplex::BlockerMapConstIterator BlockerMapConstIterator;
-//
-//	typedef typename SkeletonBlockerComplex::Visitor Visitor;
-//
-//
-//	/** @name Constructors / Destructors / Initialization
-//	 */
-//	//@{
-//
-//	explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL) :
-//    				  Skeleton_blocker_complex<SkeletonBlockerDS>(num_vertices_, visitor_) { }
-//
-//	/**
-//	 * @brief Constructor with a list of simplices
-//	 * @details The list of simplices must be the list
-//	 * of simplices of a simplicial complex, sorted with increasing dimension.
-//	 */
-//	//soon deprecated
-//	explicit Skeleton_blocker_complex(std::list<Simplex_handle>& simplices, Visitor* visitor_ = NULL) :
-//    				  Skeleton_blocker_complex<SkeletonBlockerDS>(simplices, visitor_) { }
-//
-//	/**
-//	 * @brief Constructor with a list of simplices
-//	 * @details The list of simplices must be the list of simplices of a simplicial complex.
-//	 */
-//	template<typename SimpleHandleOutputIterator>
-//	explicit Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin,SimpleHandleOutputIterator simplex_end,bool is_flag_complex = false,Visitor* visitor_ = NULL) :
-//	Skeleton_blocker_complex<SkeletonBlockerDS>(simplex_begin,simplex_end,is_flag_complex, visitor_) { }
-//
-//
-//	virtual ~Skeleton_blocker_complex() { }
-//
-//	//@}
-//
-//	/**
-//	 * Returns true iff the blocker 'sigma' is popable.
-//	 * To define popable, let us call 'L' the complex that
-//	 * consists in the current complex without the blocker 'sigma'.
-//	 * A blocker 'sigma' is then "popable" if the link of 'sigma'
-//	 * in L is reducible.
-//	 *
-//	 */
-//	bool is_popable_blocker(Blocker_handle sigma) const;
-//
-//	/**
-//	 * Removes all the popable blockers of the complex and delete them.
-//	 * @returns the number of popable blockers deleted
-//	 */
-//	void remove_popable_blockers();
-//
-//	/**
-//	 * Removes all the popable blockers of the complex passing through v and delete them.
-//	 */
-//	void remove_popable_blockers(Vertex_handle v);
-//
-//	/**
-//	 * @brief Removes all the popable blockers of the complex passing through v and delete them.
-//	 * Also remove popable blockers in the neighborhood if they became popable.
-//	 *
-//	 */
-//	void remove_all_popable_blockers(Vertex_handle v);
-//
-//	/**
-//	 * Remove the star of the vertex 'v'
-//	 */
-//	void remove_star(Vertex_handle v) ;
-//
-//private:
-//	/**
-//	 * after removing the star of a simplex, blockers sigma that contains this simplex must be removed.
-//	 * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added
-//	 * whenever the dimension of tau is at least 2.
-//	 */
-//	void update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed);
-//
-//public:
-//	/**
-//	 * Remove the star of the edge connecting vertices a and b.
-//	 * @returns the number of blocker that have been removed
-//	 */
-//	void remove_star(Vertex_handle a, Vertex_handle b) ;
-//
-//	/**
-//	 * Remove the star of the edge 'e'.
-//	 */
-//	void remove_star(Edge_handle e) ;
-//
-//	/**
-//	 * Remove the star of the simplex 'sigma' which needs to belong to the complex
-//	 */
-//	void remove_star(const Simplex_handle& sigma);
-//
-//	/**
-//	 * @brief add a maximal simplex plus all its cofaces.
-//	 * @details the simplex must have dimension greater than one (otherwise use add_vertex or add_edge).
-//	 */
-//	void add_simplex(const Simplex_handle& sigma);
-//
-//private:
-//	/**
-//	 * remove all blockers that contains sigma
-//	 */
-//	void remove_blocker_containing_simplex(const Simplex_handle& sigma) ;
-//
-//	/**
-//	 * remove all blockers that contains sigma
-//	 */
-//	void remove_blocker_include_in_simplex(const Simplex_handle& sigma);
-//
-//public:
-//	enum simplifiable_status {
-//		NOT_HOMOTOPY_EQ, MAYBE_HOMOTOPY_EQ, HOMOTOPY_EQ
-//	};
-//
-//	simplifiable_status is_remove_star_homotopy_preserving(const Simplex_handle& simplex) {
-//		// todo write a virtual method 'link' in Skeleton_blocker_complex which will be overloaded by the current one of Skeleton_blocker_geometric_complex
-//		// then call it there to build the link and return the value of link.is_contractible()
-//		return MAYBE_HOMOTOPY_EQ;
-//	}
-//
-//	enum contractible_status {
-//		NOT_CONTRACTIBLE, MAYBE_CONTRACTIBLE, CONTRACTIBLE
-//	};
-//
-//	/**
-//	 * @brief %Test if the complex is reducible using a strategy defined in the class
-//	 * (by default it tests if the complex is a cone)
-//	 * @details Note that NO could be returned if some invariant ensures that the complex
-//	 * is not a point (for instance if the euler characteristic is different from 1).
-//	 * This function will surely have to return MAYBE in some case because the
-//	 * associated problem is undecidable but it in practice, it can often
-//	 * be solved with the help of geometry.
-//	 */
-//	virtual contractible_status is_contractible() const {
-//		if (this->is_cone()) {
-//			return CONTRACTIBLE;
-//		} else {
-//			return MAYBE_CONTRACTIBLE;
-//		}
-//	}
-//
-//	/** @Edge contraction operations
-//	 */
-//	//@{
-//
-//	/**
-//	 * @return If ignore_popable_blockers is true
-//	 * then the result is true iff the link condition at edge ab is satisfied
-//	 * or equivalently iff no blocker contains ab.
-//	 * If ignore_popable_blockers is false then the
-//	 * result is true iff all blocker containing ab are popable.
-//	 */
-//	bool link_condition(Vertex_handle a, Vertex_handle b, bool ignore_popable_blockers = false) const{
-//		for (auto blocker : this->const_blocker_range(a))
-//			if (blocker->contains(b)) {
-//				// false if ignore_popable_blockers is false
-//				// otherwise the blocker has to be popable
-//				return ignore_popable_blockers && is_popable_blocker(blocker);
-//			}
-//		return true;
-//
-//	}
-//
-//	/**
-//	 * @return If ignore_popable_blockers is true
-//	 * then the result is true iff the link condition at edge ab is satisfied
-//	 * or equivalently iff no blocker contains ab.
-//	 * If ignore_popable_blockers is false then the
-//	 * result is true iff all blocker containing ab are popable.
-//	 */
-//	bool link_condition(Edge_handle e, bool ignore_popable_blockers = false) const {
-//		const Graph_edge& edge = (*this)[e];
-//		assert(this->get_address(edge.first()));
-//		assert(this->get_address(edge.second()));
-//		Vertex_handle a(*this->get_address(edge.first()));
-//		Vertex_handle b(*this->get_address(edge.second()));
-//		return link_condition(a, b, ignore_popable_blockers);
-//	}
-//
-//protected:
-//	/**
-//	 * Compute simplices beta such that a.beta is an order 0 blocker
-//	 * that may be used to construct a new blocker after contracting ab.
-//	 * It requires that the link condition is satisfied.
-//	 */
-//	void tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const;
-//
-//private:
-//	/**
-//	 * @brief "Replace" the edge 'bx' by the edge 'ax'.
-//	 * Assume that the edge 'bx' was present whereas 'ax' was not.
-//	 * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'.
-//	 * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added
-//	 * and just before edge 'bx' had been removed. That way, it can
-//	 * eventually access to information of 'ax'.
-//	 */
-//	void swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x);
-//
-//private:
-//	/**
-//	 * @brief removes all blockers passing through the edge 'ab'
-//	 */
-//	void delete_blockers_around_vertex(Vertex_handle v);
-//
-//	/**
-//	 * @brief removes all blockers passing through the edge 'ab'
-//	 */
-//	void delete_blockers_around_edge(Vertex_handle a, Vertex_handle b);
-//
-//public:
-//	/**
-//	 * Contracts the edge.
-//	 * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
-//	 * it removes first all blockers passing through 'ab'
-//	 */
-//	void contract_edge(Edge_handle edge) {
-//		contract_edge(this->first_vertex(edge), this->second_vertex(edge));
-//	}
-//
-//	/**
-//	 * Contracts the edge connecting vertices a and b.
-//	 * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
-//	 * it removes first all blockers passing through 'ab'
-//	 */
-//	void contract_edge(Vertex_handle a, Vertex_handle b);
-//
-//private:
-//	void get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, std::set<Simplex_handle>& blockers_to_add);
-//
-//	/**
-//	 * delete all blockers that passes through a or b
-//	 */
-//	void delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b);
-//	void update_edges_after_contraction(Vertex_handle a, Vertex_handle b) ;
-//
-//	void notify_changed_edges(Vertex_handle a) ;
-//	//@}
-//};
-
-
-
+/**
+ * Returns true iff the blocker 'sigma' is popable.
+ * To define popable, let us call 'L' the complex that
+ * consists in the current complex without the blocker 'sigma'.
+ * A blocker 'sigma' is then "popable" if the link of 'sigma'
+ * in L is reducible.
+ *
+ */
 template<typename SkeletonBlockerDS>
 bool Skeleton_blocker_complex<SkeletonBlockerDS>::is_popable_blocker(Blocker_handle sigma) const {
-	assert(this->contains_blocker(*sigma));
-	Skeleton_blocker_link_complex<Skeleton_blocker_complex> link_blocker_sigma;
-	build_link_of_blocker(*this, *sigma, link_blocker_sigma);
+  assert(this->contains_blocker(*sigma));
+  Skeleton_blocker_link_complex<Skeleton_blocker_complex> link_blocker_sigma;
+  build_link_of_blocker(*this, *sigma, link_blocker_sigma);
 
-	bool res = link_blocker_sigma.is_contractible() == CONTRACTIBLE;
-	return res;
+  bool res = link_blocker_sigma.is_contractible() == CONTRACTIBLE;
+  return res;
 }
 
-
 /**
  * Removes all the popable blockers of the complex and delete them.
  * @returns the number of popable blockers deleted
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers() {
-	std::list<Vertex_handle> vertex_to_check;
-	for (auto v : this->vertex_range())
-		vertex_to_check.push_front(v);
-
-	while (!vertex_to_check.empty()) {
-		Vertex_handle v = vertex_to_check.front();
-		vertex_to_check.pop_front();
-
-		bool blocker_popable_found = true;
-		while (blocker_popable_found) {
-			blocker_popable_found = false;
-			for (auto block : this->blocker_range(v)) {
-				if (this->is_popable_blocker(block)) {
-					for (Vertex_handle nv : *block)
-						if (nv != v) vertex_to_check.push_back(nv);
-					this->delete_blocker(block);
-					blocker_popable_found = true;
-					break;
-				}
-			}
-		}
-	}
+  std::list<Vertex_handle> vertex_to_check;
+  for (auto v : this->vertex_range())
+    vertex_to_check.push_front(v);
+
+  while (!vertex_to_check.empty()) {
+    Vertex_handle v = vertex_to_check.front();
+    vertex_to_check.pop_front();
+
+    bool blocker_popable_found = true;
+    while (blocker_popable_found) {
+      blocker_popable_found = false;
+      for (auto block : this->blocker_range(v)) {
+        if (this->is_popable_blocker(block)) {
+          for (Vertex_handle nv : *block)
+            if (nv != v) vertex_to_check.push_back(nv);
+          this->delete_blocker(block);
+          blocker_popable_found = true;
+          break;
+        }
+      }
+    }
+  }
 }
 
 /**
@@ -349,16 +85,16 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers() {
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers(Vertex_handle v) {
-	bool blocker_popable_found = true;
-	while (blocker_popable_found) {
-		blocker_popable_found = false;
-		for (auto block : this->blocker_range(v)) {
-			if (is_popable_blocker(block)) {
-				this->delete_blocker(block);
-				blocker_popable_found = true;
-			}
-		}
-	}
+  bool blocker_popable_found = true;
+  while (blocker_popable_found) {
+    blocker_popable_found = false;
+    for (auto block : this->blocker_range(v)) {
+      if (is_popable_blocker(block)) {
+        this->delete_blocker(block);
+        blocker_popable_found = true;
+      }
+    }
+  }
 }
 
 /**
@@ -368,73 +104,83 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers(Vertex
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_all_popable_blockers(Vertex_handle v) {
-	std::list<Vertex_handle> vertex_to_check;
-	vertex_to_check.push_front(v);
-
-	while (!vertex_to_check.empty()) {
-		Vertex_handle v = vertex_to_check.front();
-		vertex_to_check.pop_front();
-
-		bool blocker_popable_found = true;
-		while (blocker_popable_found) {
-			blocker_popable_found = false;
-			for (auto block : this->blocker_range(v)) {
-				if (this->is_popable_blocker(block)) {
-					for (Vertex_handle nv : *block)
-						if (nv != v) vertex_to_check.push_back(nv);
-					this->delete_blocker(block);
-					blocker_popable_found = true;
-					break;
-				}
-			}
-		}
-	}
+  std::list<Vertex_handle> vertex_to_check;
+  vertex_to_check.push_front(v);
+
+  while (!vertex_to_check.empty()) {
+    Vertex_handle v = vertex_to_check.front();
+    vertex_to_check.pop_front();
+
+    bool blocker_popable_found = true;
+    while (blocker_popable_found) {
+      blocker_popable_found = false;
+      for (auto block : this->blocker_range(v)) {
+        if (this->is_popable_blocker(block)) {
+          for (Vertex_handle nv : *block)
+            if (nv != v) vertex_to_check.push_back(nv);
+          this->delete_blocker(block);
+          blocker_popable_found = true;
+          break;
+        }
+      }
+    }
+  }
 }
 
-
-
+/**
+ * Remove the star of the vertice 'v'
+ */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle v) {
-	// we remove the blockers that are not consistent anymore
-	update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(v));
-	while (this->degree(v) > 0) {
-		Vertex_handle w(* (adjacent_vertices(v.vertex, this->skeleton).first));
-		this->remove_edge(v, w);
-	}
-	this->remove_vertex(v);
+  // we remove the blockers that are not consistent anymore
+  update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(v));
+  while (this->degree(v) > 0) {
+    Vertex_handle w(* (adjacent_vertices(v.vertex, this->skeleton).first));
+    this->remove_edge(v, w);
+  }
+  this->remove_vertex(v);
 }
 
+/**
+ * after removing the star of a simplex, blockers sigma that contains this simplex must be removed.
+ * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added
+ * whenever the dimension of tau is at least 2.
+ */
 template<typename SkeletonBlockerDS>
-void Skeleton_blocker_complex<SkeletonBlockerDS>::update_blockers_after_remove_star_of_vertex_or_edge(const Simplex_handle& simplex_to_be_removed) {
-	std::list <Blocker_handle> blockers_to_update;
-	if (simplex_to_be_removed.empty()) return;
-
-	auto v0 = simplex_to_be_removed.first_vertex();
-	for (auto blocker : this->blocker_range(v0)) {
-		if (blocker->contains(simplex_to_be_removed))
-			blockers_to_update.push_back(blocker);
-	}
-
-	for (auto blocker_to_update : blockers_to_update) {
-		Simplex_handle sub_blocker_to_be_added;
-		bool sub_blocker_need_to_be_added =
-				(blocker_to_update->dimension() - simplex_to_be_removed.dimension()) >= 2;
-		if (sub_blocker_need_to_be_added) {
-			sub_blocker_to_be_added = *blocker_to_update;
-			sub_blocker_to_be_added.difference(simplex_to_be_removed);
-		}
-		this->delete_blocker(blocker_to_update);
-		if (sub_blocker_need_to_be_added)
-			this->add_blocker(sub_blocker_to_be_added);
-	}
+void Skeleton_blocker_complex<SkeletonBlockerDS>::update_blockers_after_remove_star_of_vertex_or_edge(
+    const Simplex_handle& simplex_to_be_removed) {
+  std::list <Blocker_handle> blockers_to_update;
+  if (simplex_to_be_removed.empty()) return;
+
+  auto v0 = simplex_to_be_removed.first_vertex();
+  for (auto blocker : this->blocker_range(v0)) {
+    if (blocker->contains(simplex_to_be_removed))
+      blockers_to_update.push_back(blocker);
+  }
+
+  for (auto blocker_to_update : blockers_to_update) {
+    Simplex_handle sub_blocker_to_be_added;
+    bool sub_blocker_need_to_be_added =
+        (blocker_to_update->dimension() - simplex_to_be_removed.dimension()) >= 2;
+    if (sub_blocker_need_to_be_added) {
+      sub_blocker_to_be_added = *blocker_to_update;
+      sub_blocker_to_be_added.difference(simplex_to_be_removed);
+    }
+    this->delete_blocker(blocker_to_update);
+    if (sub_blocker_need_to_be_added)
+      this->add_blocker(sub_blocker_to_be_added);
+  }
 }
 
-
+/**
+ * Remove the star of the edge connecting vertices a and b.
+ * @returns the number of blocker that have been removed
+ */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle a, Vertex_handle b) {
-	update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(a, b));
-	// we remove the edge
-	this->remove_edge(a, b);
+  update_blockers_after_remove_star_of_vertex_or_edge(Simplex_handle(a, b));
+  // we remove the edge
+  this->remove_edge(a, b);
 }
 
 /**
@@ -442,22 +188,23 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle a, V
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Edge_handle e) {
-	return remove_star(this->first_vertex(e), this->second_vertex(e));
+  return remove_star(this->first_vertex(e), this->second_vertex(e));
 }
 
-
-
+/**
+ * Remove the star of the simplex 'sigma' which needs to belong to the complex
+ */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(const Simplex_handle& sigma) {
-	assert(this->contains(sigma));
-	if (sigma.dimension() == 0) {
-		remove_star(sigma.first_vertex());
-	} else if (sigma.dimension() == 1) {
-		remove_star(sigma.first_vertex(), sigma.last_vertex());
-	} else {
-		remove_blocker_containing_simplex(sigma);
-		this->add_blocker(sigma);
-	}
+  assert(this->contains(sigma));
+  if (sigma.dimension() == 0) {
+    remove_star(sigma.first_vertex());
+  } else if (sigma.dimension() == 1) {
+    remove_star(sigma.first_vertex(), sigma.last_vertex());
+  } else {
+    remove_blocker_containing_simplex(sigma);
+    this->add_blocker(sigma);
+  }
 }
 
 /**
@@ -466,35 +213,34 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(const Simplex_hand
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::add_simplex(const Simplex_handle& sigma) {
-	assert(!this->contains(sigma));
-	assert(sigma.dimension() > 1);
-
-	int num_vertex_to_add = 0;
-	for(auto v : sigma)
-		if(!contains_vertex(v)) ++num_vertex_to_add;
-	while(num_vertex_to_add--) add_vertex();
-
-	for(auto u_it = sigma.begin(); u_it != sigma.end(); ++u_it)
-		for(auto v_it = u_it; ++v_it != sigma.end(); /**/){
-			std::cout <<"add edge"<<*u_it<<" "<<*v_it<<std::endl;
-			add_edge(*u_it,*v_it);
-		}
-	remove_blocker_include_in_simplex(sigma);
+  assert(!this->contains(sigma));
+  assert(sigma.dimension() > 1);
+
+  int num_vertex_to_add = 0;
+  for (auto v : sigma)
+    if (!contains_vertex(v)) ++num_vertex_to_add;
+  while (num_vertex_to_add--) add_vertex();
+
+  for (auto u_it = sigma.begin(); u_it != sigma.end(); ++u_it)
+    for (auto v_it = u_it; ++v_it != sigma.end(); /**/) {
+      std::cout << "add edge" << *u_it << " " << *v_it << std::endl;
+      add_edge(*u_it, *v_it);
+    }
+  remove_blocker_include_in_simplex(sigma);
 }
 
-
 /**
  * remove all blockers that contains sigma
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_containing_simplex(const Simplex_handle& sigma) {
-	std::vector <Blocker_handle> blockers_to_remove;
-	for (auto blocker : this->blocker_range(sigma.first_vertex())) {
-		if (blocker->contains(sigma))
-			blockers_to_remove.push_back(blocker);
-	}
-	for (auto blocker_to_update : blockers_to_remove)
-		this->delete_blocker(blocker_to_update);
+  std::vector <Blocker_handle> blockers_to_remove;
+  for (auto blocker : this->blocker_range(sigma.first_vertex())) {
+    if (blocker->contains(sigma))
+      blockers_to_remove.push_back(blocker);
+  }
+  for (auto blocker_to_update : blockers_to_remove)
+    this->delete_blocker(blocker_to_update);
 }
 
 /**
@@ -502,173 +248,195 @@ void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_containing_simp
  */
 template<typename SkeletonBlockerDS>
 void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_include_in_simplex(const Simplex_handle& sigma) {
-	std::vector <Blocker_handle> blockers_to_remove;
-	for (auto blocker : this->blocker_range(sigma.first_vertex())) {
-		if (sigma.contains(*blocker))
-			blockers_to_remove.push_back(blocker);
-	}
-	for (auto blocker_to_update : blockers_to_remove)
-		this->delete_blocker(blocker_to_update);
+  std::vector <Blocker_handle> blockers_to_remove;
+  for (auto blocker : this->blocker_range(sigma.first_vertex())) {
+    if (sigma.contains(*blocker))
+      blockers_to_remove.push_back(blocker);
+  }
+  for (auto blocker_to_update : blockers_to_remove)
+    this->delete_blocker(blocker_to_update);
 }
 
-
+/**
+ * Compute simplices beta such that a.beta is an order 0 blocker
+ * that may be used to construct a new blocker after contracting ab.
+ * It requires that the link condition is satisfied.
+ */
 template<typename SkeletonBlockerDS>
-void Skeleton_blocker_complex<SkeletonBlockerDS>::tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex_handle> & buffer) const {
-	for (auto const & blocker : this->const_blocker_range(a)) {
-		Simplex_handle beta = (*blocker);
-		beta.remove_vertex(a);
-		buffer.push_back(beta);
-	}
-
-	Simplex_handle n;
-	this->add_neighbours(b, n);
-	this->remove_neighbours(a, n);
-	n.remove_vertex(a);
-
-
-	for (Vertex_handle y : n) {
-		Simplex_handle beta;
-		beta.add_vertex(y);
-		buffer.push_back(beta);
-	}
+void Skeleton_blocker_complex<SkeletonBlockerDS>::tip_blockers(Vertex_handle a, Vertex_handle b,
+                                                               std::vector<Simplex_handle> & buffer) const {
+  for (auto const & blocker : this->const_blocker_range(a)) {
+    Simplex_handle beta = (*blocker);
+    beta.remove_vertex(a);
+    buffer.push_back(beta);
+  }
+
+  Simplex_handle n;
+  this->add_neighbours(b, n);
+  this->remove_neighbours(a, n);
+  n.remove_vertex(a);
+
+
+  for (Vertex_handle y : n) {
+    Simplex_handle beta;
+    beta.add_vertex(y);
+    buffer.push_back(beta);
+  }
 }
 
+/**
+ * @brief "Replace" the edge 'bx' by the edge 'ax'.
+ * Assume that the edge 'bx' was present whereas 'ax' was not.
+ * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'.
+ * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added
+ * and just before edge 'bx' had been removed. That way, it can
+ * eventually access to information of 'ax'.
+ */
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x) {
-	this->add_edge(a, x);
-	if (this->visitor) this->visitor->on_swaped_edge(a, b, x);
-	this->remove_edge(b, x);
+  this->add_edge(a, x);
+  if (this->visitor) this->visitor->on_swaped_edge(a, b, x);
+  this->remove_edge(b, x);
 }
 
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_vertex(Vertex_handle v) {
-	std::list <Blocker_handle> blockers_to_delete;
-	for (auto blocker : this->blocker_range(v)) {
-		blockers_to_delete.push_back(blocker);
-	}
-	while (!blockers_to_delete.empty()) {
-		this->remove_blocker(blockers_to_delete.back());
-		blockers_to_delete.pop_back();
-	}
+  std::list <Blocker_handle> blockers_to_delete;
+  for (auto blocker : this->blocker_range(v)) {
+    blockers_to_delete.push_back(blocker);
+  }
+  while (!blockers_to_delete.empty()) {
+    this->remove_blocker(blockers_to_delete.back());
+    blockers_to_delete.pop_back();
+  }
 }
 
+/**
+ * @brief removes all blockers passing through the edge 'ab'
+ */
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_edge(Vertex_handle a, Vertex_handle b) {
-	std::list<Blocker_handle> blocker_to_delete;
-	for (auto blocker : this->blocker_range(a))
-		if (blocker->contains(b)) blocker_to_delete.push_back(blocker);
-	while (!blocker_to_delete.empty()) {
-		this->delete_blocker(blocker_to_delete.back());
-		blocker_to_delete.pop_back();
-	}
+  std::list<Blocker_handle> blocker_to_delete;
+  for (auto blocker : this->blocker_range(a))
+    if (blocker->contains(b)) blocker_to_delete.push_back(blocker);
+  while (!blocker_to_delete.empty()) {
+    this->delete_blocker(blocker_to_delete.back());
+    blocker_to_delete.pop_back();
+  }
 }
 
+/**
+ * Contracts the edge connecting vertices a and b.
+ * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied,
+ * it removes first all blockers passing through 'ab'
+ */
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::contract_edge(Vertex_handle a, Vertex_handle b) {
-	assert(this->contains_vertex(a));
-	assert(this->contains_vertex(b));
-	assert(this->contains_edge(a, b));
+  assert(this->contains_vertex(a));
+  assert(this->contains_vertex(b));
+  assert(this->contains_edge(a, b));
 
-	// if some blockers passes through 'ab', we remove them.
-	if (!link_condition(a, b))
-		delete_blockers_around_edge(a, b);
+  // if some blockers passes through 'ab', we remove them.
+  if (!link_condition(a, b))
+    delete_blockers_around_edge(a, b);
 
-	std::set<Simplex_handle> blockers_to_add;
+  std::set<Simplex_handle> blockers_to_add;
 
-	get_blockers_to_be_added_after_contraction(a, b, blockers_to_add);
+  get_blockers_to_be_added_after_contraction(a, b, blockers_to_add);
 
-	delete_blockers_around_vertices(a, b);
+  delete_blockers_around_vertices(a, b);
 
-	update_edges_after_contraction(a, b);
+  update_edges_after_contraction(a, b);
 
-	this->remove_vertex(b);
+  this->remove_vertex(b);
 
-	notify_changed_edges(a);
+  notify_changed_edges(a);
 
-	for (auto block : blockers_to_add)
-		this->add_blocker(block);
+  for (auto block : blockers_to_add)
+    this->add_blocker(block);
 
-	assert(this->contains_vertex(a));
-	assert(!this->contains_vertex(b));
+  assert(this->contains_vertex(a));
+  assert(!this->contains_vertex(b));
 }
 
 template<typename SkeletonBlockerDS>
 void
-Skeleton_blocker_complex<SkeletonBlockerDS>::get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, std::set<Simplex_handle>& blockers_to_add) {
-	blockers_to_add.clear();
-
-	typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > LinkComplexType;
-
-	LinkComplexType link_a(*this, a);
-	LinkComplexType link_b(*this, b);
-
-	std::vector<Simplex_handle> vector_alpha, vector_beta;
-
-	tip_blockers(a, b, vector_alpha);
-	tip_blockers(b, a, vector_beta);
-
-	for (auto alpha = vector_alpha.begin(); alpha != vector_alpha.end(); ++alpha) {
-		for (auto beta = vector_beta.begin(); beta != vector_beta.end(); ++beta) {
-			Simplex_handle sigma = *alpha;
-			sigma.union_vertices(*beta);
-			Root_simplex_handle sigma_id = this->get_id(sigma);
-			if (this->contains(sigma) &&
-					proper_faces_in_union<Skeleton_blocker_complex<SkeletonBlockerDS>>(sigma_id, link_a, link_b)) {
-				// Blocker_handle blocker = new Simplex_handle(sigma);
-				sigma.add_vertex(a);
-				blockers_to_add.insert(sigma);
-			}
-		}
-	}
+Skeleton_blocker_complex<SkeletonBlockerDS>::get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, 
+    std::set<Simplex_handle>& blockers_to_add) {
+  blockers_to_add.clear();
+
+  typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > LinkComplexType;
+
+  LinkComplexType link_a(*this, a);
+  LinkComplexType link_b(*this, b);
+
+  std::vector<Simplex_handle> vector_alpha, vector_beta;
+
+  tip_blockers(a, b, vector_alpha);
+  tip_blockers(b, a, vector_beta);
+
+  for (auto alpha = vector_alpha.begin(); alpha != vector_alpha.end(); ++alpha) {
+    for (auto beta = vector_beta.begin(); beta != vector_beta.end(); ++beta) {
+      Simplex_handle sigma = *alpha;
+      sigma.union_vertices(*beta);
+      Root_simplex_handle sigma_id = this->get_id(sigma);
+      if (this->contains(sigma) &&
+          proper_faces_in_union<Skeleton_blocker_complex < SkeletonBlockerDS >> (sigma_id, link_a, link_b)) {
+        // Blocker_handle blocker = new Simplex_handle(sigma);
+        sigma.add_vertex(a);
+        blockers_to_add.insert(sigma);
+      }
+    }
+  }
 }
 
+/**
+ * delete all blockers that passes through a or b
+ */
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b) {
-	std::vector<Blocker_handle> blocker_to_delete;
-	for (auto bl : this->blocker_range(a))
-		blocker_to_delete.push_back(bl);
-	for (auto bl : this->blocker_range(b))
-		blocker_to_delete.push_back(bl);
-	while (!blocker_to_delete.empty()) {
-		this->delete_blocker(blocker_to_delete.back());
-		blocker_to_delete.pop_back();
-	}
+  std::vector<Blocker_handle> blocker_to_delete;
+  for (auto bl : this->blocker_range(a))
+    blocker_to_delete.push_back(bl);
+  for (auto bl : this->blocker_range(b))
+    blocker_to_delete.push_back(bl);
+  while (!blocker_to_delete.empty()) {
+    this->delete_blocker(blocker_to_delete.back());
+    blocker_to_delete.pop_back();
+  }
 }
 
-
-
 template<typename SkeletonBlockerDS>
 void
 Skeleton_blocker_complex<SkeletonBlockerDS>::update_edges_after_contraction(Vertex_handle a, Vertex_handle b) {
-	// We update the set of edges
-	this->remove_edge(a, b);
-
-	// For all edges {b,x} incident to b,
-	// we remove {b,x} and add {a,x} if not already there.
-	while (this->degree(b) > 0) {
-		Vertex_handle x(*(adjacent_vertices(b.vertex, this->skeleton).first));
-		if (!this->contains_edge(a, x))
-			// we 'replace' the edge 'bx' by the edge 'ax'
-			this->swap_edge(a, b, x);
-		else
-			this->remove_edge(b, x);
-	}
+  // We update the set of edges
+  this->remove_edge(a, b);
+
+  // For all edges {b,x} incident to b,
+  // we remove {b,x} and add {a,x} if not already there.
+  while (this->degree(b) > 0) {
+    Vertex_handle x(*(adjacent_vertices(b.vertex, this->skeleton).first));
+    if (!this->contains_edge(a, x))
+      // we 'replace' the edge 'bx' by the edge 'ax'
+      this->swap_edge(a, b, x);
+    else
+      this->remove_edge(b, x);
+  }
 }
 
 template<typename SkeletonBlockerDS>
 void
-Skeleton_blocker_complex<SkeletonBlockerDS>::notify_changed_edges(Vertex_handle a){
-	// We notify the visitor that all edges incident to 'a' had changed
-	boost_adjacency_iterator v, v_end;
-
-	for (tie(v, v_end) = adjacent_vertices(a.vertex, this->skeleton); v != v_end; ++v)
-		if (this->visitor) this->visitor->on_changed_edge(a, Vertex_handle(*v));
+Skeleton_blocker_complex<SkeletonBlockerDS>::notify_changed_edges(Vertex_handle a) {
+  // We notify the visitor that all edges incident to 'a' had changed
+  boost_adjacency_iterator v, v_end;
 
+  for (tie(v, v_end) = adjacent_vertices(a.vertex, this->skeleton); v != v_end; ++v)
+    if (this->visitor) this->visitor->on_changed_edge(a, Vertex_handle(*v));
 }
 
 
diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h
index 1a96128f..c6968e44 100644
--- a/src/Witness_complex/include/gudhi/Witness_complex.h
+++ b/src/Witness_complex/include/gudhi/Witness_complex.h
@@ -273,16 +273,12 @@ private:
     if (!map.empty())
       {
         std::cout << map.begin()->first;
-        if (map.begin()->second.children() == root())
-          std::cout << "Sweet potato";
         if (has_children(map.begin()))
           print_sc(map.begin()->second.children());
         typename Dictionary::iterator it;
         for (it = map.begin()+1; it != map.end(); ++it)
           {
             std::cout << "," << it->first;
-            if (map.begin()->second.children() == root())
-              std::cout << "Sweet potato";
             if (has_children(it))
               print_sc(it->second.children());
           }
diff --git a/src/Witness_complex/include/gudhi/Witness_complex1.h b/src/Witness_complex/include/gudhi/Witness_complex1.h
index 6950c3bf..49ba7529 100644
--- a/src/Witness_complex/include/gudhi/Witness_complex1.h
+++ b/src/Witness_complex/include/gudhi/Witness_complex1.h
@@ -334,10 +334,10 @@ private:
                 if (j != i)
                   {
                     std::cout << "+++ We are at vertex=" << knn[witness_id][j] << std::endl;
-                    if (curr_sibl->find(knn[witness_id][j]) == null_simplex())
+                    if (curr_sibl->members().find(knn[witness_id][j]) == null_simplex())
                       return false;
                     std::cout << "++++ the simplex is there\n";
-                    curr_sh = curr_sibl->find(knn[witness_id][j]);
+                    curr_sh = curr_sibl->members().find(knn[witness_id][j]);
                     std::cout << "++++ curr_sh affectation is OK\n";
                     if (has_children(curr_sh))
                       curr_sibl = curr_sh->second.children();