integrated kernel code in pers_representation branch

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

Former-commit-id: 58e77263a0da3674e4699cef832b6d357dcf12e2

7 files changed, 246 insertions, 207 deletions

diff --git a/src/Simplex_tree/doc/Intro_simplex_tree.h b/src/Simplex_tree/doc/Intro_simplex_tree.h
index 769491d9..6b80d1c9 100644
--- a/src/Simplex_tree/doc/Intro_simplex_tree.h
+++ b/src/Simplex_tree/doc/Intro_simplex_tree.h
@@ -79,7 +79,6 @@ Number of vertices = 10   Number of simplices = 98 \endcode
  * 1 incidence relations in a complex. It is consequently faster when accessing the boundary of a simplex, but is less
  * compact and harder to construct from scratch.
  * 
- * \copyright GNU General Public License v3.
  * @}
  */
 
diff --git a/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp b/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp
index 217e251f..9bd51106 100644
--- a/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp
+++ b/src/Simplex_tree/example/cech_complex_cgal_mini_sphere_3d.cpp
@@ -1,5 +1,5 @@
-/*    This file is part of the Gudhi Library. The Gudhi library 
- *    (Geometric Understanding in Higher Dimensions) is a generic C++ 
+/*    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):       Clément Maria
@@ -33,7 +33,7 @@
 
 #include <string>
 #include <vector>
-#include <limits>  // infinity
+#include <limits>   // infinity
 #include <utility>  // for pair
 #include <map>
 
@@ -50,15 +50,14 @@ using Vertex_handle = Simplex_tree::Vertex_handle;
 using Simplex_handle = Simplex_tree::Simplex_handle;
 using Filtration_value = Simplex_tree::Filtration_value;
 using Siblings = Simplex_tree::Siblings;
-using Graph_t = boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS
-, boost::property < Gudhi::vertex_filtration_t, Filtration_value >
-, boost::property < Gudhi::edge_filtration_t, Filtration_value >
->;
-using Edge_t = std::pair< Vertex_handle, Vertex_handle >;
+using Graph_t = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS,
+                                      boost::property<Gudhi::vertex_filtration_t, Filtration_value>,
+                                      boost::property<Gudhi::edge_filtration_t, Filtration_value> >;
+using Edge_t = std::pair<Vertex_handle, Vertex_handle>;
 
-using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<3> >;
+using Kernel = CGAL::Epick_d<CGAL::Dimension_tag<3> >;
 using Point = Kernel::Point_d;
-using Traits = CGAL::Min_sphere_of_points_d_traits_d<Kernel,Filtration_value,3>;
+using Traits = CGAL::Min_sphere_of_points_d_traits_d<Kernel, Filtration_value, 3>;
 using Min_sphere = CGAL::Min_sphere_of_spheres_d<Traits>;
 
 using Points_off_reader = Gudhi::Points_off_reader<Point>;
@@ -76,7 +75,7 @@ class Cech_blocker {
       std::cout << vertex << ", ";
 #endif  // DEBUG_TRACES
     }
-    Min_sphere ms(points.begin(),points.end());
+    Min_sphere ms(points.begin(), points.end());
     Filtration_value radius = ms.radius();
 #if DEBUG_TRACES
     std::cout << "] - radius = " << radius << " - returns " << (radius > threshold_) << std::endl;
@@ -85,24 +84,20 @@ class Cech_blocker {
     return (radius > threshold_);
   }
   Cech_blocker(Simplex_tree& simplex_tree, Filtration_value threshold, const std::vector<Point>& point_cloud)
-    : simplex_tree_(simplex_tree),
-      threshold_(threshold),
-      point_cloud_(point_cloud) { }
+      : simplex_tree_(simplex_tree), threshold_(threshold), point_cloud_(point_cloud) {}
+
  private:
   Simplex_tree simplex_tree_;
   Filtration_value threshold_;
   std::vector<Point> point_cloud_;
 };
 
-template< typename InputPointRange>
-Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold);
+template <typename InputPointRange>
+Graph_t compute_proximity_graph(InputPointRange& points, Filtration_value threshold);
 
-void program_options(int argc, char * argv[]
-                     , std::string & off_file_points
-                     , Filtration_value & threshold
-                     , int & dim_max);
+void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& threshold, int& dim_max);
 
-int main(int argc, char * argv[]) {
+int main(int argc, char* argv[]) {
   std::string off_file_points;
   Filtration_value threshold;
   int dim_max;
@@ -115,7 +110,7 @@ int main(int argc, char * argv[]) {
   // Compute the proximity graph of the points
   Graph_t prox_graph = compute_proximity_graph(off_reader.get_point_cloud(), threshold);
 
-  //Min_sphere sph1(off_reader.get_point_cloud()[0], off_reader.get_point_cloud()[1], off_reader.get_point_cloud()[2]);
+  // Min_sphere sph1(off_reader.get_point_cloud()[0], off_reader.get_point_cloud()[1], off_reader.get_point_cloud()[2]);
   // Construct the Rips complex in a Simplex Tree
   Simplex_tree st;
   // insert the proximity graph in the simplex tree
@@ -135,7 +130,8 @@ int main(int argc, char * argv[]) {
   std::cout << "* The complex contains " << st.num_simplices() << " simplices - dimension=" << st.dimension() << "\n";
   std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
   for (auto f_simplex : st.filtration_simplex_range()) {
-    std::cout << "   " << "[" << st.filtration(f_simplex) << "] ";
+    std::cout << "   "
+              << "[" << st.filtration(f_simplex) << "] ";
     for (auto vertex : st.simplex_vertex_range(f_simplex)) {
       std::cout << static_cast<int>(vertex) << " ";
     }
@@ -145,24 +141,19 @@ int main(int argc, char * argv[]) {
   return 0;
 }
 
-void program_options(int argc, char * argv[]
-                     , std::string & off_file_points
-                     , Filtration_value & threshold
-                     , int & dim_max) {
+void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& threshold, int& dim_max) {
   namespace po = boost::program_options;
   po::options_description hidden("Hidden options");
-  hidden.add_options()
-      ("input-file", po::value<std::string>(&off_file_points),
-       "Name of an OFF file containing a 3d point set.\n");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of an OFF file containing a 3d point set.\n");
 
   po::options_description visible("Allowed options", 100);
-  visible.add_options()
-      ("help,h", "produce help message")
-      ("max-edge-length,r",
-       po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()),
-       "Maximal length of an edge for the Cech complex construction.")
-      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
-       "Maximal dimension of the Cech complex we want to compute.");
+  visible.add_options()("help,h", "produce help message")(
+      "max-edge-length,r",
+      po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()),
+      "Maximal length of an edge for the Cech complex construction.")(
+      "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
+      "Maximal dimension of the Cech complex we want to compute.");
 
   po::positional_options_description pos;
   pos.add("input-file", 1);
@@ -171,8 +162,7 @@ void program_options(int argc, char * argv[]
   all.add(visible).add(hidden);
 
   po::variables_map vm;
-  po::store(po::command_line_parser(argc, argv).
-            options(all).positional(pos).run(), vm);
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
   po::notify(vm);
 
   if (vm.count("help") || !vm.count("input-file")) {
@@ -194,10 +184,10 @@ void program_options(int argc, char * argv[]
  * The type PointCloud furnishes .begin() and .end() methods, that return
  * iterators with value_type Point.
  */
-template< typename InputPointRange>
-Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold) {
-  std::vector< Edge_t > edges;
-  std::vector< Filtration_value > edges_fil;
+template <typename InputPointRange>
+Graph_t compute_proximity_graph(InputPointRange& points, Filtration_value threshold) {
+  std::vector<Edge_t> edges;
+  std::vector<Filtration_value> edges_fil;
 
   Kernel k;
   Vertex_handle idx_u, idx_v;
@@ -217,16 +207,13 @@ Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value thresh
     ++idx_u;
   }
 
-  Graph_t skel_graph(edges.begin()
-                     , edges.end()
-                     , edges_fil.begin()
-                     , idx_u);  // number of points labeled from 0 to idx_u-1
+  Graph_t skel_graph(edges.begin(), edges.end(), edges_fil.begin(),
+                     idx_u);  // number of points labeled from 0 to idx_u-1
 
   auto vertex_prop = boost::get(Gudhi::vertex_filtration_t(), skel_graph);
 
   boost::graph_traits<Graph_t>::vertex_iterator vi, vi_end;
-  for (std::tie(vi, vi_end) = boost::vertices(skel_graph);
-       vi != vi_end; ++vi) {
+  for (std::tie(vi, vi_end) = boost::vertices(skel_graph); vi != vi_end; ++vi) {
     boost::put(vertex_prop, *vi, 0.);
   }
 
diff --git a/src/Simplex_tree/example/graph_expansion_with_blocker.cpp b/src/Simplex_tree/example/graph_expansion_with_blocker.cpp
index 86bfb8cb..0d458cbd 100644
--- a/src/Simplex_tree/example/graph_expansion_with_blocker.cpp
+++ b/src/Simplex_tree/example/graph_expansion_with_blocker.cpp
@@ -27,8 +27,7 @@
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Simplex_handle = Simplex_tree::Simplex_handle;
 
-int main(int argc, char * const argv[]) {
-
+int main(int argc, char* const argv[]) {
   // Construct the Simplex Tree with a 1-skeleton graph example
   Simplex_tree simplexTree;
 
@@ -45,33 +44,32 @@ int main(int argc, char * const argv[]) {
   simplexTree.insert_simplex({5, 6}, 10.);
   simplexTree.insert_simplex({6}, 10.);
 
-  simplexTree.expansion_with_blockers(3, [&](Simplex_handle sh){
-      bool result = false;
-      std::cout << "Blocker on [";
-      // User can loop on the vertices from the given simplex_handle i.e.
-      for (auto vertex : simplexTree.simplex_vertex_range(sh)) {
-        // We block the expansion, if the vertex '6' is in the given list of vertices
-        if (vertex == 6)
-          result = true;
-        std::cout << vertex << ", ";
-      }
-      std::cout << "] ( " << simplexTree.filtration(sh);
-      // User can re-assign a new filtration value directly in the blocker (default is the maximal value of boudaries)
-      simplexTree.assign_filtration(sh, simplexTree.filtration(sh) + 1.);
+  simplexTree.expansion_with_blockers(3, [&](Simplex_handle sh) {
+    bool result = false;
+    std::cout << "Blocker on [";
+    // User can loop on the vertices from the given simplex_handle i.e.
+    for (auto vertex : simplexTree.simplex_vertex_range(sh)) {
+      // We block the expansion, if the vertex '6' is in the given list of vertices
+      if (vertex == 6) result = true;
+      std::cout << vertex << ", ";
+    }
+    std::cout << "] ( " << simplexTree.filtration(sh);
+    // User can re-assign a new filtration value directly in the blocker (default is the maximal value of boudaries)
+    simplexTree.assign_filtration(sh, simplexTree.filtration(sh) + 1.);
 
-      std::cout << " + 1. ) = " << result << std::endl;
+    std::cout << " + 1. ) = " << result << std::endl;
 
-      return result;
-    });
+    return result;
+  });
 
   std::cout << "********************************************************************\n";
   std::cout << "* The complex contains " << simplexTree.num_simplices() << " simplices";
   std::cout << "   - dimension " << simplexTree.dimension() << "\n";
   std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
   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 << "(" << vertex << ")";
+    std::cout << "   "
+              << "[" << simplexTree.filtration(f_simplex) << "] ";
+    for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) std::cout << "(" << vertex << ")";
     std::cout << std::endl;
   }
 
diff --git a/src/Simplex_tree/example/simple_simplex_tree.cpp b/src/Simplex_tree/example/simple_simplex_tree.cpp
index b6b65b88..828977c2 100644
--- a/src/Simplex_tree/example/simple_simplex_tree.cpp
+++ b/src/Simplex_tree/example/simple_simplex_tree.cpp
@@ -30,10 +30,10 @@
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Vertex_handle = Simplex_tree::Vertex_handle;
 using Filtration_value = Simplex_tree::Filtration_value;
-using typeVectorVertex = std::vector< Vertex_handle >;
-using typePairSimplexBool = std::pair< Simplex_tree::Simplex_handle, bool >;
+using typeVectorVertex = std::vector<Vertex_handle>;
+using typePairSimplexBool = std::pair<Simplex_tree::Simplex_handle, bool>;
 
-int main(int argc, char * const argv[]) {
+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;
@@ -54,7 +54,7 @@ int main(int argc, char * const argv[]) {
 
   // ++ FIRST
   std::cout << "   * INSERT 0" << std::endl;
-  typeVectorVertex firstSimplexVector = { 0 };
+  typeVectorVertex firstSimplexVector = {0};
   typePairSimplexBool returnValue =
       simplexTree.insert_simplex(firstSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
 
@@ -66,9 +66,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ SECOND
   std::cout << "   * INSERT 1" << std::endl;
-  typeVectorVertex secondSimplexVector = { 1 };
-  returnValue =
-      simplexTree.insert_simplex(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
+  typeVectorVertex secondSimplexVector = {1};
+  returnValue = simplexTree.insert_simplex(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + 1 INSERTED" << std::endl;
@@ -78,9 +77,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ THIRD
   std::cout << "   * INSERT (0,1)" << std::endl;
-  typeVectorVertex thirdSimplexVector = { 0, 1 };
-  returnValue =
-      simplexTree.insert_simplex(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  typeVectorVertex thirdSimplexVector = {0, 1};
+  returnValue = simplexTree.insert_simplex(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (0,1) INSERTED" << std::endl;
@@ -90,9 +88,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ FOURTH
   std::cout << "   * INSERT 2" << std::endl;
-  typeVectorVertex fourthSimplexVector = { 2 };
-  returnValue =
-      simplexTree.insert_simplex(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
+  typeVectorVertex fourthSimplexVector = {2};
+  returnValue = simplexTree.insert_simplex(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + 2 INSERTED" << std::endl;
@@ -102,9 +99,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ FIFTH
   std::cout << "   * INSERT (2,0)" << std::endl;
-  typeVectorVertex fifthSimplexVector = { 2, 0 };
-  returnValue =
-      simplexTree.insert_simplex(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  typeVectorVertex fifthSimplexVector = {2, 0};
+  returnValue = simplexTree.insert_simplex(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (2,0) INSERTED" << std::endl;
@@ -114,9 +110,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ SIXTH
   std::cout << "   * INSERT (2,1)" << std::endl;
-  typeVectorVertex sixthSimplexVector = { 2, 1 };
-  returnValue =
-      simplexTree.insert_simplex(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  typeVectorVertex sixthSimplexVector = {2, 1};
+  returnValue = simplexTree.insert_simplex(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (2,1) INSERTED" << std::endl;
@@ -126,9 +121,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ SEVENTH
   std::cout << "   * INSERT (2,1,0)" << std::endl;
-  typeVectorVertex seventhSimplexVector = { 2, 1, 0 };
-  returnValue =
-      simplexTree.insert_simplex(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE));
+  typeVectorVertex seventhSimplexVector = {2, 1, 0};
+  returnValue = simplexTree.insert_simplex(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (2,1,0) INSERTED" << std::endl;
@@ -138,9 +132,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ EIGHTH
   std::cout << "   * INSERT 3" << std::endl;
-  typeVectorVertex eighthSimplexVector = { 3 };
-  returnValue =
-      simplexTree.insert_simplex(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
+  typeVectorVertex eighthSimplexVector = {3};
+  returnValue = simplexTree.insert_simplex(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + 3 INSERTED" << std::endl;
@@ -150,9 +143,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ NINETH
   std::cout << "   * INSERT (3,0)" << std::endl;
-  typeVectorVertex ninethSimplexVector = { 3, 0 };
-  returnValue =
-      simplexTree.insert_simplex(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  typeVectorVertex ninethSimplexVector = {3, 0};
+  returnValue = simplexTree.insert_simplex(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (3,0) INSERTED" << std::endl;
@@ -162,7 +154,7 @@ int main(int argc, char * const argv[]) {
 
   // ++ TENTH
   std::cout << "   * INSERT 0 (already inserted)" << std::endl;
-  typeVectorVertex tenthSimplexVector = { 0 };
+  typeVectorVertex tenthSimplexVector = {0};
   // With a different filtration value
   returnValue = simplexTree.insert_simplex(tenthSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE));
 
@@ -174,9 +166,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ ELEVENTH
   std::cout << "   * INSERT (2,1,0) (already inserted)" << std::endl;
-  typeVectorVertex eleventhSimplexVector = { 2, 1, 0 };
-  returnValue =
-      simplexTree.insert_simplex(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE));
+  typeVectorVertex eleventhSimplexVector = {2, 1, 0};
+  returnValue = simplexTree.insert_simplex(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (2,1,0) INSERTED" << std::endl;
@@ -192,9 +183,9 @@ int main(int argc, char * const argv[]) {
   std::cout << "   - dimension " << simplexTree.dimension() << "\n";
   std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
   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 << "(" << vertex << ")";
+    std::cout << "   "
+              << "[" << simplexTree.filtration(f_simplex) << "] ";
+    for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) std::cout << "(" << vertex << ")";
     std::cout << std::endl;
   }
   //   [0.1] 0
@@ -217,7 +208,7 @@ int main(int argc, char * const argv[]) {
   else
     std::cout << "***- NO IT ISN'T\n";
 
-  typeVectorVertex unknownSimplexVector = { 15 };
+  typeVectorVertex unknownSimplexVector = {15};
   simplexFound = simplexTree.find(unknownSimplexVector);
   std::cout << "**************IS THE SIMPLEX {15} IN THE SIMPLEX TREE ?\n";
   if (simplexFound != simplexTree.null_simplex())
@@ -232,7 +223,7 @@ int main(int argc, char * const argv[]) {
   else
     std::cout << "***- NO IT ISN'T\n";
 
-  typeVectorVertex otherSimplexVector = { 1, 15 };
+  typeVectorVertex otherSimplexVector = {1, 15};
   simplexFound = simplexTree.find(otherSimplexVector);
   std::cout << "**************IS THE SIMPLEX {15,1} IN THE SIMPLEX TREE ?\n";
   if (simplexFound != simplexTree.null_simplex())
@@ -240,7 +231,7 @@ int main(int argc, char * const argv[]) {
   else
     std::cout << "***- NO IT ISN'T\n";
 
-  typeVectorVertex invSimplexVector = { 1, 2, 0 };
+  typeVectorVertex invSimplexVector = {1, 2, 0};
   simplexFound = simplexTree.find(invSimplexVector);
   std::cout << "**************IS THE SIMPLEX {1,2,0} IN THE SIMPLEX TREE ?\n";
   if (simplexFound != simplexTree.null_simplex())
@@ -248,7 +239,7 @@ int main(int argc, char * const argv[]) {
   else
     std::cout << "***- NO IT ISN'T\n";
 
-  simplexFound = simplexTree.find({ 0, 1 });
+  simplexFound = simplexTree.find({0, 1});
   std::cout << "**************IS THE SIMPLEX {0,1} IN THE SIMPLEX TREE ?\n";
   if (simplexFound != simplexTree.null_simplex())
     std::cout << "***+ YES IT IS!\n";
@@ -256,23 +247,20 @@ int main(int argc, char * const argv[]) {
     std::cout << "***- NO IT ISN'T\n";
 
   std::cout << "**************COFACES OF {0,1} IN CODIMENSION 1 ARE\n";
-  for (auto& simplex : simplexTree.cofaces_simplex_range(simplexTree.find({0,1}), 1)) {
-    for (auto vertex : simplexTree.simplex_vertex_range(simplex))
-      std::cout << "(" << vertex << ")";
+  for (auto& simplex : simplexTree.cofaces_simplex_range(simplexTree.find({0, 1}), 1)) {
+    for (auto vertex : simplexTree.simplex_vertex_range(simplex)) std::cout << "(" << vertex << ")";
     std::cout << std::endl;
   }
 
   std::cout << "**************STARS OF {0,1} ARE\n";
-  for (auto& simplex : simplexTree.star_simplex_range(simplexTree.find({0,1}))) {
-    for (auto vertex : simplexTree.simplex_vertex_range(simplex))
-      std::cout << "(" << vertex << ")";
+  for (auto& simplex : simplexTree.star_simplex_range(simplexTree.find({0, 1}))) {
+    for (auto vertex : simplexTree.simplex_vertex_range(simplex)) std::cout << "(" << vertex << ")";
     std::cout << std::endl;
   }
 
   std::cout << "**************BOUNDARIES OF {0,1,2} ARE\n";
-  for (auto& simplex : simplexTree.boundary_simplex_range(simplexTree.find({0,1,2}))) {
-    for (auto vertex : simplexTree.simplex_vertex_range(simplex))
-      std::cout << "(" << vertex << ")";
+  for (auto& simplex : simplexTree.boundary_simplex_range(simplexTree.find({0, 1, 2}))) {
+    for (auto vertex : simplexTree.simplex_vertex_range(simplex)) std::cout << "(" << vertex << ")";
     std::cout << std::endl;
   }
 
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index cb6ab309..7456cb1f 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -49,6 +49,7 @@
 #include <initializer_list>
 #include <algorithm>  // for std::max
 #include <cstdint>  // for std::uint32_t
+#include <iterator>  // for std::distance
 
 namespace Gudhi {
 
@@ -106,8 +107,9 @@ class Simplex_tree {
   };
   struct Key_simplex_base_dummy {
     Key_simplex_base_dummy() {}
-    void assign_key(Simplex_key) { }
-    Simplex_key key() const { assert(false); return -1; }
+    // Undefined so it will not link
+    void assign_key(Simplex_key);
+    Simplex_key key() const;
   };
   typedef typename std::conditional<Options::store_key, Key_simplex_base_real, Key_simplex_base_dummy>::type
       Key_simplex_base;
@@ -121,7 +123,7 @@ class Simplex_tree {
   };
   struct Filtration_simplex_base_dummy {
     Filtration_simplex_base_dummy() {}
-    void assign_filtration(Filtration_value f) { assert(f == 0); }
+    void assign_filtration(Filtration_value GUDHI_CHECK_code(f)) { GUDHI_CHECK(f == 0, "filtration value specified for a complex that does not store them"); }
     Filtration_value filtration() const { return 0; }
   };
   typedef typename std::conditional<Options::store_filtration, Filtration_simplex_base_real,
@@ -500,6 +502,7 @@ class Simplex_tree {
    * sh has children.*/
   template<class SimplexHandle>
   bool has_children(SimplexHandle sh) const {
+    // Here we rely on the root using null_vertex(), which cannot match any real vertex.
     return (sh->second.children()->parent() == sh->first);
   }
 
@@ -529,18 +532,30 @@ class Simplex_tree {
   Simplex_handle find_simplex(const std::vector<Vertex_handle> & simplex) {
     Siblings * tmp_sib = &root_;
     Dictionary_it tmp_dit;
-    Vertex_handle last = simplex.back();
-    for (auto v : simplex) {
-      tmp_dit = tmp_sib->members_.find(v);
-      if (tmp_dit == tmp_sib->members_.end()) {
+    auto vi = simplex.begin();
+    if (Options::contiguous_vertices) {
+      // Equivalent to the first iteration of the normal loop
+      GUDHI_CHECK(contiguous_vertices(), "non-contiguous vertices");
+      Vertex_handle v = *vi++;
+      if(v < 0 || v >= static_cast<Vertex_handle>(root_.members_.size()))
         return null_simplex();
-      }
-      if (!has_children(tmp_dit) && v != last) {
+      tmp_dit = root_.members_.begin() + v;
+      if (vi == simplex.end())
+        return tmp_dit;
+      if (!has_children(tmp_dit))
+        return null_simplex();
+      tmp_sib = tmp_dit->second.children();
+    }
+    for (;;) {
+      tmp_dit = tmp_sib->members_.find(*vi++);
+      if (tmp_dit == tmp_sib->members_.end())
+        return null_simplex();
+      if (vi == simplex.end())
+        return tmp_dit;
+      if (!has_children(tmp_dit))
         return null_simplex();
-      }
       tmp_sib = tmp_dit->second.children();
     }
-    return tmp_dit;
   }
 
   /** \brief Returns the Simplex_handle corresponding to the 0-simplex
@@ -584,12 +599,14 @@ class Simplex_tree {
     std::pair<Simplex_handle, bool> res_insert;
     auto vi = simplex.begin();
     for (; vi != simplex.end() - 1; ++vi) {
+      GUDHI_CHECK(*vi != null_vertex(), "cannot use the dummy null_vertex() as a real vertex");
       res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration));
       if (!(has_children(res_insert.first))) {
         res_insert.first->second.assign_children(new Siblings(curr_sib, *vi));
       }
       curr_sib = res_insert.first->second.children();
     }
+    GUDHI_CHECK(*vi != null_vertex(), "cannot use the dummy null_vertex() as a real vertex");
     res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration));
     if (!res_insert.second) {
       // if already in the complex
@@ -664,71 +681,67 @@ class Simplex_tree {
    */
   template<class InputVertexRange = std::initializer_list<Vertex_handle>>
   std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(const InputVertexRange& Nsimplex,
-                                   Filtration_value filtration = 0) {
+				    Filtration_value filtration = 0) {
     auto first = std::begin(Nsimplex);
     auto last = std::end(Nsimplex);
 
     if (first == last)
-      return std::pair<Simplex_handle, bool>(null_simplex(), true);  // ----->>
+      return { null_simplex(), true };  // ----->>
 
     // Copy before sorting
-    std::vector<Vertex_handle> copy(first, last);
+    thread_local std::vector<Vertex_handle> copy;
+    copy.clear();
+    copy.insert(copy.end(), first, last);
     std::sort(std::begin(copy), std::end(copy));
+    GUDHI_CHECK_code(
+      for (Vertex_handle v : copy)
+        GUDHI_CHECK(v != null_vertex(), "cannot use the dummy null_vertex() as a real vertex");
+    )
 
-    std::vector<std::vector<Vertex_handle>> to_be_inserted;
-    std::vector<std::vector<Vertex_handle>> to_be_propagated;
-    return rec_insert_simplex_and_subfaces(copy, to_be_inserted, to_be_propagated, filtration);
+    return insert_simplex_and_subfaces_sorted(copy, filtration);
   }
 
  private:
-  std::pair<Simplex_handle, bool> rec_insert_simplex_and_subfaces(std::vector<Vertex_handle>& the_simplex,
-                                                             std::vector<std::vector<Vertex_handle>>& to_be_inserted,
-                                                             std::vector<std::vector<Vertex_handle>>& to_be_propagated,
-                                                             Filtration_value filtration = 0.0) {
-    std::pair<Simplex_handle, bool> insert_result;
-    if (the_simplex.size() > 1) {
-      // Get and remove last vertex
-      Vertex_handle last_vertex = the_simplex.back();
-      the_simplex.pop_back();
-      // Recursive call after last vertex removal
-      insert_result = rec_insert_simplex_and_subfaces(the_simplex, to_be_inserted, to_be_propagated, filtration);
-
-      // Concatenation of to_be_inserted and to_be_propagated
-      to_be_inserted.insert(to_be_inserted.begin(), to_be_propagated.begin(), to_be_propagated.end());
-      to_be_propagated = to_be_inserted;
-
-      // to_be_inserted treatment
-      for (auto& simplex_tbi : to_be_inserted) {
-        simplex_tbi.push_back(last_vertex);
-      }
-      std::vector<Vertex_handle> last_simplex(1, last_vertex);
-      to_be_inserted.insert(to_be_inserted.begin(), last_simplex);
-      // i.e. (0,1,2) =>
-      // [to_be_inserted | to_be_propagated] = [(1) (0,1) | (0)]
-      // [to_be_inserted | to_be_propagated] = [(2) (0,2) (1,2) (0,1,2) | (0) (1) (0,1)]
-      // N.B. : it is important the last inserted to be the highest in dimension
-      // in order to return the "last" insert_simplex result
-
-      // insert all to_be_inserted
-      for (auto& simplex_tbi : to_be_inserted) {
-        insert_result = insert_vertex_vector(simplex_tbi, filtration);
-      }
-    } else if (the_simplex.size() == 1) {
-      // When reaching the end of recursivity, vector of simplices shall be empty and filled on back recursive
-      if ((to_be_inserted.size() != 0) || (to_be_propagated.size() != 0)) {
-        std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Error vector not empty\n";
-        exit(-1);
+  /// Same as insert_simplex_and_subfaces but assumes that the range of vertices is sorted
+  template<class ForwardVertexRange = std::initializer_list<Vertex_handle>>
+  std::pair<Simplex_handle, bool> insert_simplex_and_subfaces_sorted(const ForwardVertexRange& Nsimplex, Filtration_value filt = 0) {
+    auto first = std::begin(Nsimplex);
+    auto last = std::end(Nsimplex);
+    if (first == last)
+      return { null_simplex(), true }; // FIXME: false would make more sense to me.
+    GUDHI_CHECK(std::is_sorted(first, last), "simplex vertices listed in unsorted order");
+    // Update dimension if needed. We could wait to see if the insertion succeeds, but I doubt there is much to gain.
+    dimension_ = (std::max)(dimension_, static_cast<int>(std::distance(first, last)) - 1);
+    return rec_insert_simplex_and_subfaces_sorted(root(), first, last, filt);
+  }
+  // To insert {1,2,3,4}, we insert {2,3,4} twice, once at the root, and once below 1.
+  template<class ForwardVertexIterator>
+  std::pair<Simplex_handle, bool> rec_insert_simplex_and_subfaces_sorted(Siblings* sib, ForwardVertexIterator first, ForwardVertexIterator last, Filtration_value filt) {
+    // An alternative strategy would be:
+    // - try to find the complete simplex, if found (and low filtration) exit
+    // - insert all the vertices at once in sib
+    // - loop over those (new or not) simplices, with a recursive call(++first, last)
+    Vertex_handle vertex_one = *first;
+    auto&& dict = sib->members();
+    auto insertion_result = dict.emplace(vertex_one, Node(sib, filt));
+    Simplex_handle simplex_one = insertion_result.first;
+    bool one_is_new = insertion_result.second;
+    if (!one_is_new) {
+      if (filtration(simplex_one) > filt) {
+	assign_filtration(simplex_one, filt);
+      } else {
+	// FIXME: this interface makes no sense, and it doesn't seem to be tested.
+	insertion_result.first = null_simplex();
       }
-      std::vector<Vertex_handle> first_simplex(1, the_simplex.back());
-      // i.e. (0,1,2) => [to_be_inserted | to_be_propagated] = [(0) | ]
-      to_be_inserted.push_back(first_simplex);
-
-      insert_result = insert_vertex_vector(first_simplex, filtration);
-    } else {
-        std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Recursivity error\n";
-        exit(-1);
     }
-    return insert_result;
+    if (++first == last) return insertion_result;
+    if (!has_children(simplex_one))
+      // TODO: have special code here, we know we are building the whole subtree from scratch.
+      simplex_one->second.assign_children(new Siblings(sib, vertex_one));
+    auto res = rec_insert_simplex_and_subfaces_sorted(simplex_one->second.children(), first, last, filt);
+    // No need to continue if the full simplex was already there with a low enough filtration value.
+    if (res.first != null_simplex()) rec_insert_simplex_and_subfaces_sorted(sib, first, last, filt);
+    return res;
   }
 
  public:
@@ -941,8 +954,9 @@ class Simplex_tree {
    * called.
    *
    * Inserts all vertices and edges given by a OneSkeletonGraph.
-   * OneSkeletonGraph must be a model of boost::AdjacencyGraph,
-   * boost::EdgeListGraph and boost::PropertyGraph.
+   * OneSkeletonGraph must be a model of
+   * <a href="http://www.boost.org/doc/libs/1_65_1/libs/graph/doc/EdgeListGraph.html">boost::EdgeListGraph</a>
+   * and <a href="http://www.boost.org/doc/libs/1_65_1/libs/graph/doc/PropertyGraph.html">boost::PropertyGraph</a>.
    *
    * The vertex filtration value is accessible through the property tag
    * vertex_filtration_t.
@@ -952,7 +966,10 @@ class Simplex_tree {
    * boost::graph_traits<OneSkeletonGraph>::vertex_descriptor
    *                                    must be Vertex_handle.
    * boost::graph_traits<OneSkeletonGraph>::directed_category
-   *                                    must be undirected_tag. */
+   *                                    must be undirected_tag.
+   *
+   * If an edge appears with multiplicity, the function will arbitrarily pick
+   * one representative to read the filtration value.  */
   template<class OneSkeletonGraph>
   void insert_graph(const OneSkeletonGraph& skel_graph) {
     // the simplex tree must be empty
@@ -983,18 +1000,22 @@ class Simplex_tree {
          ++e_it) {
       auto u = source(*e_it, skel_graph);
       auto v = target(*e_it, skel_graph);
-      if (u < v) {
-        // count edges only once { std::swap(u,v); } // u < v
-        auto sh = find_vertex(u);
-        if (!has_children(sh)) {
-          sh->second.assign_children(new Siblings(&root_, sh->first));
-        }
-
-        sh->second.children()->members().emplace(
-                                                 v,
-                                                 Node(sh->second.children(),
-                                                      boost::get(edge_filtration_t(), skel_graph, *e_it)));
+      if (u == v) throw "Self-loops are not simplicial";
+      // We cannot skip edges with the wrong orientation and expect them to
+      // come a second time with the right orientation, that does not always
+      // happen in practice. emplace() should be a NOP when an element with the
+      // same key is already there, so seeing the same edge multiple times is
+      // ok.
+      // Should we actually forbid multiple edges? That would be consistent
+      // with rejecting self-loops.
+      if (v < u) std::swap(u, v);
+      auto sh = find_vertex(u);
+      if (!has_children(sh)) {
+        sh->second.assign_children(new Siblings(&root_, sh->first));
       }
+
+      sh->second.children()->members().emplace(v,
+          Node(sh->second.children(), boost::get(edge_filtration_t(), skel_graph, *e_it)));
     }
   }
 
@@ -1138,7 +1159,7 @@ class Simplex_tree {
             to_be_inserted=false;
             break;
           }
-          filt = std::max(filt, filtration(border_child));
+          filt = (std::max)(filt, filtration(border_child));
         }
         if (to_be_inserted) {
           intersection.emplace_back(next->first, Node(nullptr, filt));
@@ -1334,7 +1355,7 @@ class Simplex_tree {
       if (sh_dimension >= dimension_)
         // Stop browsing as soon as the dimension is reached, no need to go furter
         return false;
-      new_dimension = std::max(new_dimension, sh_dimension);
+      new_dimension = (std::max)(new_dimension, sh_dimension);
     }
     dimension_ = new_dimension;
     return true;
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h
index 7e0a454d..ab7346d4 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_iterators.h
@@ -23,6 +23,8 @@
 #ifndef SIMPLEX_TREE_SIMPLEX_TREE_ITERATORS_H_
 #define SIMPLEX_TREE_SIMPLEX_TREE_ITERATORS_H_
 
+#include <gudhi/Debug_utils.h>
+
 #include <boost/iterator/iterator_facade.hpp>
 #include <boost/version.hpp>
 #if BOOST_VERSION >= 105600
@@ -109,11 +111,18 @@ class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade<
       : last_(sh->first),
         sib_(nullptr),
         st_(st) {
+    // Only check once at the beginning instead of for every increment, as this is expensive.
+    if (SimplexTree::Options::contiguous_vertices)
+      GUDHI_CHECK(st_->contiguous_vertices(), "The set of vertices is not { 0, ..., n } without holes");
     Siblings * sib = st->self_siblings(sh);
     next_ = sib->parent();
     sib_ = sib->oncles();
     if (sib_ != nullptr) {
-      sh_ = sib_->find(next_);
+      if (SimplexTree::Options::contiguous_vertices && sib_->oncles() == nullptr)
+        // Only relevant for edges
+        sh_ = sib_->members_.begin()+next_;
+      else
+        sh_ = sib_->find(next_);
     } else {
       sh_ = st->null_simplex();
     }  // vertex: == end()
@@ -140,14 +149,19 @@ class Simplex_tree_boundary_simplex_iterator : public boost::iterator_facade<
     Siblings * for_sib = sib_;
     Siblings * new_sib = sib_->oncles();
     auto rit = suffix_.rbegin();
-    if (SimplexTree::Options::contiguous_vertices && new_sib == nullptr && rit != suffix_.rend()) {
-      // We reached the root, use a short-cut to find a vertex. We could also
-      // optimize finding the second vertex of a segment, but people are
-      // expected to call endpoints().
-      assert(st_->contiguous_vertices());
-      sh_ = for_sib->members_.begin()+*rit;
-      for_sib = sh_->second.children();
-      ++rit;
+    if (SimplexTree::Options::contiguous_vertices && new_sib == nullptr) {
+      // We reached the root, use a short-cut to find a vertex.
+      if (rit == suffix_.rend()) {
+        // Segment, this vertex is the last boundary simplex
+        sh_ = for_sib->members_.begin()+last_;
+        sib_ = nullptr;
+        return;
+      } else {
+        // Dim >= 2, initial step of the descent
+        sh_ = for_sib->members_.begin()+*rit;
+        for_sib = sh_->second.children();
+        ++rit;
+      }
     }
     for (; rit != suffix_.rend(); ++rit) {
       sh_ = for_sib->find(*rit);
diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
index 580d610a..f63ea080 100644
--- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp
+++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
@@ -863,3 +863,35 @@ BOOST_AUTO_TEST_CASE(make_filtration_non_decreasing) {
   BOOST_CHECK(st == st_other_copy);
   
 }
+
+BOOST_AUTO_TEST_CASE(insert_graph) {
+  std::cout << "********************************************************************" << std::endl;
+  std::cout << "INSERT GRAPH" << std::endl;
+  typedef typename boost::adjacency_list<boost::vecS, boost::vecS,
+          boost::undirectedS,
+          boost::property<vertex_filtration_t, double>,
+          boost::property<edge_filtration_t, double>> Graph;
+  Graph g(3);
+  // filtration value 0 everywhere
+  put(Gudhi::vertex_filtration_t(), g, 0, 0);
+  put(Gudhi::vertex_filtration_t(), g, 1, 0);
+  put(Gudhi::vertex_filtration_t(), g, 2, 0);
+  // vertices don't always occur in sorted order
+  add_edge(0, 1, 0, g);
+  add_edge(2, 1, 0, g);
+  add_edge(2, 0, 0, g);
+
+  typedef Simplex_tree<> typeST;
+  typeST st1;
+  st1.insert_graph(g);
+  BOOST_CHECK(st1.num_simplices() == 6);
+
+  // edges can have multiplicity in the graph unless we replace the first vecS with (hash_)setS
+  add_edge(1, 0, 0, g);
+  add_edge(1, 2, 0, g);
+  add_edge(0, 2, 0, g);
+  add_edge(0, 2, 0, g);
+  typeST st2;
+  st2.insert_graph(g);
+  BOOST_CHECK(st2.num_simplices() == 6);
+}