Move bottleneck_chrono in benchmark

Add test in cmake for basic example
Move CGAL gudhi patches for bottleneck in dedicated directory
Fix cpplint syntax issue


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

Former-commit-id: 24671facf791de93dc6fd94bb39ca7362bb22959

16 files changed, 590 insertions, 2348 deletions

diff --git a/src/Bottleneck_distance/benchmark/CMakeLists.txt b/src/Bottleneck_distance/benchmark/CMakeLists.txt
new file mode 100644
index 00000000..f70dd8ff
--- /dev/null
+++ b/src/Bottleneck_distance/benchmark/CMakeLists.txt
@@ -0,0 +1,13 @@
+cmake_minimum_required(VERSION 2.6)
+project(Bottleneck_distance_benchmark)
+
+
+# requires CGAL 4.8
+# cmake -DCGAL_DIR=~/workspace/CGAL-4.8 ../../..
+if(CGAL_FOUND)
+  if (NOT CGAL_VERSION VERSION_LESS 4.8.0)
+    if (EIGEN3_FOUND)
+      add_executable ( bottleneck_chrono bottleneck_chrono.cpp )
+    endif()
+  endif ()
+endif()
diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp
index a30d42b5..456c570b 100644
--- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp
+++ b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp
@@ -4,7 +4,7 @@
  *
  *    Author:       Francois Godi
  *
- *    Copyright (C) 2015  INRIA (France)
+ *    Copyright (C) 2015  INRIA
  *
  *    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
@@ -28,35 +28,35 @@
 using namespace Gudhi::persistence_diagram;
 
 
-double upper_bound = 400.; // any real >0
+double upper_bound = 400.;  // any real > 0
 
-int main(){
-    std::ofstream objetfichier;
-    objetfichier.open("results.csv", std::ios::out);
+int main() {
+  std::ofstream result_file;
+  result_file.open("results.csv", std::ios::out);
 
-    for(int n = 1000; n<=10000; n+=1000){
-    std::uniform_real_distribution<double> unif1(0.,upper_bound);
-    std::uniform_real_distribution<double> unif2(upper_bound/1000.,upper_bound/100.);
+  for (int n = 1000; n <= 10000; n += 1000) {
+    std::uniform_real_distribution<double> unif1(0., upper_bound);
+    std::uniform_real_distribution<double> unif2(upper_bound / 1000., upper_bound / 100.);
     std::default_random_engine re;
     std::vector< std::pair<double, double> > v1, v2;
     for (int i = 0; i < n; i++) {
-        double a = unif1(re);
-        double b = unif1(re);
-        double x = unif2(re);
-        double y = unif2(re);
-        v1.emplace_back(std::min(a,b), std::max(a,b));
-        v2.emplace_back(std::min(a,b)+std::min(x,y), std::max(a,b)+std::max(x,y));
-        if(i%5==0)
-            v1.emplace_back(std::min(a,b),std::min(a,b)+x);
-        if(i%3==0)
-            v2.emplace_back(std::max(a,b),std::max(a,b)+y);
+      double a = unif1(re);
+      double b = unif1(re);
+      double x = unif2(re);
+      double y = unif2(re);
+      v1.emplace_back(std::min(a, b), std::max(a, b));
+      v2.emplace_back(std::min(a, b) + std::min(x, y), std::max(a, b) + std::max(x, y));
+      if (i % 5 == 0)
+        v1.emplace_back(std::min(a, b), std::min(a, b) + x);
+      if (i % 3 == 0)
+        v2.emplace_back(std::max(a, b), std::max(a, b) + y);
     }
     std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
     double b = bottleneck_distance(v1, v2);
     std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
-    typedef std::chrono::duration<int,std::milli> millisecs_t;
-    millisecs_t duration(std::chrono::duration_cast<millisecs_t>(end-start));
-    objetfichier << n << ";" << duration.count() << ";" << b << std::endl;
-    }
-    objetfichier.close();
+    typedef std::chrono::duration<int, std::milli> millisecs_t;
+    millisecs_t duration(std::chrono::duration_cast<millisecs_t>(end - start));
+    result_file << n << ";" << duration.count() << ";" << b << std::endl;
+  }
+  result_file.close();
 }
diff --git a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h
index ebe1123b..21187f9c 100644
--- a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h
+++ b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h
@@ -4,7 +4,7 @@
  *
  *    Author:       François Godi
  *
- *    Copyright (C) 2015  INRIA (France)
+ *    Copyright (C) 2015  INRIA
  *
  *    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
@@ -20,8 +20,8 @@
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#ifndef DOC_BOTTLENECK_DISTANCE_H_
-#define DOC_BOTTLENECK_DISTANCE_H_
+#ifndef DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_
+#define DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_
 
 // needs namespace for Doxygen to link on classes
 namespace Gudhi {
@@ -48,3 +48,4 @@ namespace bottleneck_distance {
 
 }  // namespace Gudhi
 
+#endif  // DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_
diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt
index cd53ccfc..c66623e9 100644
--- a/src/Bottleneck_distance/example/CMakeLists.txt
+++ b/src/Bottleneck_distance/example/CMakeLists.txt
@@ -8,6 +8,8 @@ if(CGAL_FOUND)
     if (EIGEN3_FOUND)
       add_executable (bottleneck_read_file_example bottleneck_read_file_example.cpp)
       add_executable (bottleneck_basic_example bottleneck_basic_example.cpp)
+
+      add_test(bottleneck_basic_example ${CMAKE_CURRENT_BINARY_DIR}/bottleneck_basic_example)
     endif()
   endif ()
 endif()
diff --git a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp
index 78e00e57..91a7302f 100644
--- a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp
+++ b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp
@@ -4,7 +4,7 @@
  *
  *    Authors:       Francois Godi, small modifications by Pawel Dlotko
  *
- *    Copyright (C) 2015  INRIA (France)
+ *    Copyright (C) 2015  INRIA
  *
  *    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
@@ -25,24 +25,24 @@
 
 int main() {
 
-    std::vector< std::pair<double, double> > v1, v2;
+  std::vector< std::pair<double, double> > v1, v2;
 
-    v1.emplace_back(2.7, 3.7);
-    v1.emplace_back(9.6, 14.);
-    v1.emplace_back(34.2, 34.974);
-    v1.emplace_back(3., std::numeric_limits<double>::infinity());
+  v1.emplace_back(2.7, 3.7);
+  v1.emplace_back(9.6, 14.);
+  v1.emplace_back(34.2, 34.974);
+  v1.emplace_back(3., std::numeric_limits<double>::infinity());
 
-    v2.emplace_back(2.8, 4.45);
-    v2.emplace_back(9.5, 14.1);
-    v2.emplace_back(3.2, std::numeric_limits<double>::infinity());
+  v2.emplace_back(2.8, 4.45);
+  v2.emplace_back(9.5, 14.1);
+  v2.emplace_back(3.2, std::numeric_limits<double>::infinity());
 
 
-    double b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2);
+  double b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2);
 
-    std::cout << "Bottleneck distance = " << b << std::endl;
+  std::cout << "Bottleneck distance = " << b << std::endl;
 
-    b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1);
+  b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1);
 
-    std::cout << "Approx bottleneck distance = " << b << std::endl;
+  std::cout << "Approx bottleneck distance = " << b << std::endl;
 
 }
diff --git a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp
index ceedccc5..4c74b66e 100644
--- a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp
+++ b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp
@@ -4,7 +4,7 @@
  *
  *    Authors:       Francois Godi, small modifications by Pawel Dlotko
  *
- *    Copyright (C) 2015  INRIA (France)
+ *    Copyright (C) 2015  INRIA
  *
  *    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
@@ -28,49 +28,42 @@
 #include <sstream>
 #include <string>
 
-std::vector< std::pair<double, double> > read_diagram_from_file( const char* filename )
-{
-    std::ifstream in;
-    in.open( filename );
-    std::vector< std::pair<double, double> > result;
-    if ( !in.is_open() )
-    {
-        std::cerr << "File : " << filename << " do not exist. The program will now terminate \n";
-        throw "File do not exist \n";
-    }
+std::vector< std::pair<double, double> > read_diagram_from_file(const char* filename) {
+  std::ifstream in;
+  in.open(filename);
+  std::vector< std::pair<double, double> > result;
+  if (!in.is_open()) {
+    std::cerr << "File : " << filename << " do not exist. The program will now terminate \n";
+    throw "File do not exist \n";
+  }
 
-    std::string line;
-    while (!in.eof())
-    {
-        getline(in,line);
-        if ( line.length() != 0 )
-        {
-            std::stringstream lineSS;
-            lineSS << line;
-            double beginn, endd;
-            lineSS >> beginn;
-            lineSS >> endd;
-            result.push_back( std::make_pair( beginn , endd ) );
-        }
+  std::string line;
+  while (!in.eof()) {
+    getline(in, line);
+    if (line.length() != 0) {
+      std::stringstream lineSS;
+      lineSS << line;
+      double beginn, endd;
+      lineSS >> beginn;
+      lineSS >> endd;
+      result.push_back(std::make_pair(beginn, endd));
     }
-    in.close();
-    return result;
-}  //read_diagram_from_file
+  }
+  in.close();
+  return result;
+}  // read_diagram_from_file
 
-int main( int argc , char** argv )
-{
-  if ( argc < 3 )
-  {
-      std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file " <<
-          "should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck" <<
-          " distance (set by default to zero). The program will now terminate \n";
+int main(int argc, char** argv) {
+  if (argc < 3) {
+    std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file " <<
+        "should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck" <<
+        " distance (set by default to zero). The program will now terminate \n";
   }
-  std::vector< std::pair< double , double > > diag1 = read_diagram_from_file( argv[1] );
-  std::vector< std::pair< double , double > > diag2 = read_diagram_from_file( argv[2] );
+  std::vector< std::pair< double, double > > diag1 = read_diagram_from_file(argv[1]);
+  std::vector< std::pair< double, double > > diag2 = read_diagram_from_file(argv[2]);
   double tolerance = 0.;
-  if ( argc == 4 )
-  {
-      tolerance = atof( argv[3] );
+  if (argc == 4) {
+    tolerance = atof(argv[3]);
   }
   double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance);
   std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl;
diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h
index 42a0d444..2b7e4767 100644
--- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h
+++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h
@@ -4,7 +4,7 @@
  *
  *    Author:       Francois Godi
  *
- *    Copyright (C) 2015  INRIA (France)
+ *    Copyright (C) 2015  INRIA
  *
  *    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
@@ -31,62 +31,65 @@ namespace Gudhi {
 namespace persistence_diagram {
 
 double bottleneck_distance_approx(Persistence_graph& g, double e) {
-    double b_lower_bound = 0.;
-    double b_upper_bound = g.diameter_bound();
-    const double alpha = std::pow(g.size(), 1./5.);
-    Graph_matching m(g);
-    Graph_matching biggest_unperfect(g);
-    while (b_upper_bound - b_lower_bound > 2*e) {
-        double step = b_lower_bound + (b_upper_bound - b_lower_bound)/alpha;
-        if(step <= b_lower_bound || step >= b_upper_bound) //Avoid precision problem
-            break;
-        m.set_r(step);
-        while (m.multi_augment()); //compute a maximum matching (in the graph corresponding to the current r)
-        if (m.perfect()) {
-            m = biggest_unperfect;
-            b_upper_bound = step;
-        } else {
-            biggest_unperfect = m;
-            b_lower_bound = step;
-        }
+  double b_lower_bound = 0.;
+  double b_upper_bound = g.diameter_bound();
+  const double alpha = std::pow(g.size(), 1. / 5.);
+  Graph_matching m(g);
+  Graph_matching biggest_unperfect(g);
+  while (b_upper_bound - b_lower_bound > 2 * e) {
+    double step = b_lower_bound + (b_upper_bound - b_lower_bound) / alpha;
+    if (step <= b_lower_bound || step >= b_upper_bound)  // Avoid precision problem
+      break;
+    m.set_r(step);
+    while (m.multi_augment());  // compute a maximum matching (in the graph corresponding to the current r)
+    if (m.perfect()) {
+      m = biggest_unperfect;
+      b_upper_bound = step;
+    } else {
+      biggest_unperfect = m;
+      b_lower_bound = step;
     }
-    return (b_lower_bound + b_upper_bound)/2.;
+  }
+  return (b_lower_bound + b_upper_bound) / 2.;
 }
 
 double bottleneck_distance_exact(Persistence_graph& g) {
-    std::vector<double> sd  = g.sorted_distances();
-    long lower_bound_i = 0;
-    long upper_bound_i = sd.size()-1;
-    const double alpha = std::pow(g.size(), 1./5.);
-    Graph_matching m(g);
-    Graph_matching biggest_unperfect(g);
-    while (lower_bound_i != upper_bound_i) {
-        long step = lower_bound_i + static_cast<long>((upper_bound_i - lower_bound_i - 1)/alpha);
-        m.set_r(sd.at(step));
-        while (m.multi_augment()); //compute a maximum matching (in the graph corresponding to the current r)
-        if (m.perfect()) {
-            m = biggest_unperfect;
-            upper_bound_i = step;
-        } else {
-            biggest_unperfect = m;
-            lower_bound_i = step + 1;
-        }
+  std::vector<double> sd = g.sorted_distances();
+  long lower_bound_i = 0;
+  long upper_bound_i = sd.size() - 1;
+  const double alpha = std::pow(g.size(), 1. / 5.);
+  Graph_matching m(g);
+  Graph_matching biggest_unperfect(g);
+  while (lower_bound_i != upper_bound_i) {
+    long step = lower_bound_i + static_cast<long> ((upper_bound_i - lower_bound_i - 1) / alpha);
+    m.set_r(sd.at(step));
+    while (m.multi_augment());  // compute a maximum matching (in the graph corresponding to the current r)
+    if (m.perfect()) {
+      m = biggest_unperfect;
+      upper_bound_i = step;
+    } else {
+      biggest_unperfect = m;
+      lower_bound_i = step + 1;
     }
-    return sd.at(lower_bound_i);
+  }
+  return sd.at(lower_bound_i);
 }
 
 /** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams (see concepts).
- *  If the last parameter e is not 0, you get an additive e-approximation, which is a lot faster to compute whatever is e.
- *  Thus, by default, e is a very small positive double, actually the smallest double possible such that the floating-point inaccuracies don't lead to a failure of the algorithm.
+ * If the last parameter e is not 0, you get an additive e-approximation, which is a lot faster to compute whatever is
+ * e.
+ * Thus, by default, e is a very small positive double, actually the smallest double possible such that the
+ * floating-point inaccuracies don't lead to a failure of the algorithm.
  *
  * \ingroup bottleneck_distance
  */
 template<typename Persistence_diagram1, typename Persistence_diagram2>
-double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=std::numeric_limits<double>::min()) {
-    Persistence_graph g(diag1, diag2, e);
-    if(g.bottleneck_alive() == std::numeric_limits<double>::infinity())
-        return std::numeric_limits<double>::infinity();
-    return std::max(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e));
+double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2,
+                           double e = std::numeric_limits<double>::min()) {
+  Persistence_graph g(diag1, diag2, e);
+  if (g.bottleneck_alive() == std::numeric_limits<double>::infinity())
+    return std::numeric_limits<double>::infinity();
+  return std::max(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e));
 }
 
 }  // namespace persistence_diagram
diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h
deleted file mode 100644
index f085b0da..00000000
--- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h
+++ /dev/null
@@ -1,582 +0,0 @@
-// Copyright (c) 2002,2011,2014 Utrecht University (The Netherlands), Max-Planck-Institute Saarbruecken (Germany).
-// All rights reserved.
-//
-// This file is part of CGAL (www.cgal.org).
-// 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.
-//
-// Licensees holding a valid commercial license may use this file in
-// accordance with the commercial license agreement provided with the software.
-//
-// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
-// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
-//
-// $URL$
-// $Id$
-//
-// Author(s)     : Hans Tangelder (<hanst@cs.uu.nl>),
-//               : Waqar Khan <wkhan@mpi-inf.mpg.de>
-
-#ifndef CGAL_KD_TREE_H
-#define CGAL_KD_TREE_H
-
-#include "Kd_tree_node.h"
-
-#include <CGAL/basic.h>
-#include <CGAL/assertions.h>
-#include <vector>
-
-#include <CGAL/algorithm.h>
-#include <CGAL/internal/Get_dimension_tag.h>
-#include <CGAL/Search_traits.h>
-
-
-#include <deque>
-#include <boost/container/deque.hpp>
-#include <boost/optional.hpp>
-
-#ifdef CGAL_HAS_THREADS
-#include <CGAL/mutex.h>
-#endif
-
-namespace CGAL {
-
-//template <class SearchTraits, class Splitter_=Median_of_rectangle<SearchTraits>, class UseExtendedNode = Tag_true >
-template <class SearchTraits, class Splitter_=Sliding_midpoint<SearchTraits>, class UseExtendedNode = Tag_true >
-class Kd_tree {
-
-public:
-  typedef SearchTraits Traits;
-  typedef Splitter_ Splitter;
-  typedef typename SearchTraits::Point_d Point_d;
-  typedef typename Splitter::Container Point_container;
-
-  typedef typename SearchTraits::FT FT;
-  typedef Kd_tree_node<SearchTraits, Splitter, UseExtendedNode > Node;
-  typedef Kd_tree_leaf_node<SearchTraits, Splitter, UseExtendedNode > Leaf_node;
-  typedef Kd_tree_internal_node<SearchTraits, Splitter, UseExtendedNode > Internal_node;
-  typedef Kd_tree<SearchTraits, Splitter> Tree;
-  typedef Kd_tree<SearchTraits, Splitter,UseExtendedNode> Self;
-
-  typedef Node* Node_handle;
-  typedef const Node* Node_const_handle;
-  typedef Leaf_node* Leaf_node_handle;
-  typedef const Leaf_node* Leaf_node_const_handle;
-  typedef Internal_node* Internal_node_handle;
-  typedef const Internal_node* Internal_node_const_handle;
-  typedef typename std::vector<const Point_d*>::const_iterator Point_d_iterator;
-  typedef typename std::vector<const Point_d*>::const_iterator Point_d_const_iterator;
-  typedef typename Splitter::Separator Separator;
-  typedef typename std::vector<Point_d>::const_iterator iterator;
-  typedef typename std::vector<Point_d>::const_iterator const_iterator;
-
-  typedef typename std::vector<Point_d>::size_type size_type;
-
-  typedef typename internal::Get_dimension_tag<SearchTraits>::Dimension D;
-
-private:
-  SearchTraits traits_;
-  Splitter split;
-
-
-  // wokaround for https://svn.boost.org/trac/boost/ticket/9332
-#if   (_MSC_VER == 1800) && (BOOST_VERSION == 105500)
-  std::deque<Internal_node> internal_nodes;
-  std::deque<Leaf_node> leaf_nodes;
-#else
-  boost::container::deque<Internal_node> internal_nodes;
-  boost::container::deque<Leaf_node> leaf_nodes;
-#endif
-
-  Node_handle tree_root;
-
-  Kd_tree_rectangle<FT,D>* bbox;
-  std::vector<Point_d> pts;
-
-  // Instead of storing the points in arrays in the Kd_tree_node
-  // we put all the data in a vector in the Kd_tree.
-  // and we only store an iterator range in the Kd_tree_node.
-  //
-  std::vector<const Point_d*> data;
-
-
-  #ifdef CGAL_HAS_THREADS
-  mutable CGAL_MUTEX building_mutex;//mutex used to protect const calls inducing build()
-  #endif
-  bool built_;
-  bool removed_;
-
-  // protected copy constructor
-  Kd_tree(const Tree& tree)
-    : traits_(tree.traits_),built_(tree.built_)
-  {};
-
-
-  // Instead of the recursive construction of the tree in the class Kd_tree_node
-  // we do this in the tree class. The advantage is that we then can optimize
-  // the allocation of the nodes.
-
-  // The leaf node
-  Node_handle
-  create_leaf_node(Point_container& c)
-  {
-    Leaf_node node(true , static_cast<unsigned int>(c.size()));
-    std::ptrdiff_t tmp = c.begin() - data.begin();
-    node.data = pts.begin() + tmp;
-
-    leaf_nodes.push_back(node);
-    Leaf_node_handle nh = &leaf_nodes.back();
-
-
-    return nh;
-  }
-
-
-  // The internal node
-
-  Node_handle
-  create_internal_node(Point_container& c, const Tag_true&)
-  {
-    return create_internal_node_use_extension(c);
-  }
-
-  Node_handle
-  create_internal_node(Point_container& c, const Tag_false&)
-  {
-    return create_internal_node(c);
-  }
-
-
-
-  // TODO: Similiar to the leaf_init function above, a part of the code should be
-  //       moved to a the class Kd_tree_node.
-  //       It is not proper yet, but the goal was to see if there is
-  //       a potential performance gain through the Compact_container
-  Node_handle
-  create_internal_node_use_extension(Point_container& c)
-  {
-    Internal_node node(false);
-    internal_nodes.push_back(node);
-    Internal_node_handle nh = &internal_nodes.back();
-
-    Separator sep;
-    Point_container c_low(c.dimension(),traits_);
-    split(sep, c, c_low);
-    nh->set_separator(sep);
-
-    int cd  = nh->cutting_dimension();
-    if(!c_low.empty()){
-      nh->lower_low_val = c_low.tight_bounding_box().min_coord(cd);
-      nh->lower_high_val = c_low.tight_bounding_box().max_coord(cd);
-    }
-    else{
-      nh->lower_low_val = nh->cutting_value();
-      nh->lower_high_val = nh->cutting_value();
-    }
-    if(!c.empty()){
-      nh->upper_low_val = c.tight_bounding_box().min_coord(cd);
-      nh->upper_high_val = c.tight_bounding_box().max_coord(cd);
-    }
-    else{
-      nh->upper_low_val = nh->cutting_value();
-      nh->upper_high_val = nh->cutting_value();
-    }
-
-    CGAL_assertion(nh->cutting_value() >= nh->lower_low_val);
-    CGAL_assertion(nh->cutting_value() <= nh->upper_high_val);
-
-    if (c_low.size() > split.bucket_size()){
-      nh->lower_ch = create_internal_node_use_extension(c_low);
-    }else{
-      nh->lower_ch = create_leaf_node(c_low);
-    }
-    if (c.size() > split.bucket_size()){
-      nh->upper_ch = create_internal_node_use_extension(c);
-    }else{
-      nh->upper_ch = create_leaf_node(c);
-    }
-
-
-
-
-    return nh;
-  }
-
-
-  // Note also that I duplicated the code to get rid if the if's for
-  // the boolean use_extension which was constant over the construction
-  Node_handle
-  create_internal_node(Point_container& c)
-  {
-    Internal_node node(false);
-    internal_nodes.push_back(node);
-    Internal_node_handle nh = &internal_nodes.back();
-    Separator sep;
-
-    Point_container c_low(c.dimension(),traits_);
-    split(sep, c, c_low);
-    nh->set_separator(sep);
-
-    if (c_low.size() > split.bucket_size()){
-      nh->lower_ch = create_internal_node(c_low);
-    }else{
-      nh->lower_ch = create_leaf_node(c_low);
-    }
-    if (c.size() > split.bucket_size()){
-      nh->upper_ch = create_internal_node(c);
-    }else{
-      nh->upper_ch = create_leaf_node(c);
-    }
-
-
-
-    return nh;
-  }
-
-
-
-public:
-
-  Kd_tree(Splitter s = Splitter(),const SearchTraits traits=SearchTraits())
-    : traits_(traits),split(s), built_(false), removed_(false)
-  {}
-
-  template <class InputIterator>
-  Kd_tree(InputIterator first, InputIterator beyond,
-          Splitter s = Splitter(),const SearchTraits traits=SearchTraits())
-    : traits_(traits),split(s), built_(false), removed_(false)
-  {
-    pts.insert(pts.end(), first, beyond);
-  }
-
-  bool empty() const {
-    return pts.empty();
-  }
-
-  void
-  build()
-  {
-    // This function is not ready to be called when a tree already exists, one
-    // must call invalidate_built() first.
-    CGAL_assertion(!is_built());
-    CGAL_assertion(!removed_);
-    const Point_d& p = *pts.begin();
-    typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits_.construct_cartesian_const_iterator_d_object();
-    int dim = static_cast<int>(std::distance(ccci(p), ccci(p,0)));
-
-    data.reserve(pts.size());
-    for(unsigned int i = 0; i < pts.size(); i++){
-      data.push_back(&pts[i]);
-    }
-    Point_container c(dim, data.begin(), data.end(),traits_);
-    bbox = new Kd_tree_rectangle<FT,D>(c.bounding_box());
-    if (c.size() <= split.bucket_size()){
-      tree_root = create_leaf_node(c);
-    }else {
-      tree_root = create_internal_node(c, UseExtendedNode());
-    }
-
-    //Reorder vector for spatial locality
-    std::vector<Point_d> ptstmp;
-    ptstmp.resize(pts.size());
-    for (std::size_t i = 0; i < pts.size(); ++i){
-      ptstmp[i] = *data[i];
-    }
-    for(std::size_t i = 0; i < leaf_nodes.size(); ++i){
-      std::ptrdiff_t tmp = leaf_nodes[i].begin() - pts.begin();
-      leaf_nodes[i].data = ptstmp.begin() + tmp;
-    }
-    pts.swap(ptstmp);
-
-    data.clear();
-
-    built_ = true;
-  }
-
-private:
-  //any call to this function is for the moment not threadsafe
-  void const_build() const {
-    #ifdef CGAL_HAS_THREADS
-    //this ensure that build() will be called once
-    CGAL_SCOPED_LOCK(building_mutex);
-    if(!is_built())
-    #endif
-      const_cast<Self*>(this)->build(); //THIS IS NOT THREADSAFE
-  }
-public:
-
-  bool is_built() const
-  {
-    return built_;
-  }
-
-  void invalidate_built()
-  {
-    if(removed_){
-      // Walk the tree to collect the remaining points.
-      // Writing directly to pts would likely work, but better be safe.
-      std::vector<Point_d> ptstmp;
-      //ptstmp.resize(root()->num_items());
-      root()->tree_items(std::back_inserter(ptstmp));
-      pts.swap(ptstmp);
-      removed_=false;
-      CGAL_assertion(is_built()); // the rest of the cleanup must happen
-    }
-    if(is_built()){
-      internal_nodes.clear();
-      leaf_nodes.clear();
-      data.clear();
-      delete bbox;
-      built_ = false;
-    }
-  }
-
-  void clear()
-  {
-    invalidate_built();
-    pts.clear();
-    removed_ = false;
-  }
-
-  void
-  insert(const Point_d& p)
-  {
-    invalidate_built();
-    pts.push_back(p);
-  }
-
-  template <class InputIterator>
-  void
-  insert(InputIterator first, InputIterator beyond)
-  {
-    invalidate_built();
-    pts.insert(pts.end(),first, beyond);
-  }
-
-private:
-  struct Equal_by_coordinates {
-    SearchTraits const* traits;
-    Point_d const* pp;
-    bool operator()(Point_d const&q) const {
-      typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits->construct_cartesian_const_iterator_d_object();
-      return std::equal(ccci(*pp), ccci(*pp,0), ccci(q));
-    }
-  };
-  Equal_by_coordinates equal_by_coordinates(Point_d const&p){
-    Equal_by_coordinates ret = { &traits(), &p };
-    return ret;
-  }
-
-public:
-  void
-  remove(const Point_d& p)
-  {
-    remove(p, equal_by_coordinates(p));
-  }
-
-  template<class Equal>
-  void
-  remove(const Point_d& p, Equal const& equal_to_p)
-  {
-#if 0
-    // This code could have quadratic runtime.
-    if (!is_built()) {
-      std::vector<Point_d>::iterator pi = std::find(pts.begin(), pts.end(), p);
-      // Precondition: the point must be there.
-      CGAL_assertion (pi != pts.end());
-      pts.erase(pi);
-      return;
-    }
-#endif
-    bool success = remove_(p, 0, false, 0, false, root(), equal_to_p);
-    CGAL_assertion(success);
-
-    // Do not set the flag is the tree has been cleared.
-    if(is_built())
-      removed_ |= success;
-  }
-private:
-  template<class Equal>
-  bool remove_(const Point_d& p,
-      Internal_node_handle grandparent, bool parent_islower,
-      Internal_node_handle parent, bool islower,
-      Node_handle node, Equal const& equal_to_p) {
-    // Recurse to locate the point
-    if (!node->is_leaf()) {
-      Internal_node_handle newparent = static_cast<Internal_node_handle>(node);
-      // FIXME: This should be if(x<y) remove low; else remove up;
-      if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] <= newparent->cutting_value()) {
-        if (remove_(p, parent, islower, newparent, true, newparent->lower(), equal_to_p))
-          return true;
-      }
-      //if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] >= newparent->cutting_value())
-        return remove_(p, parent, islower, newparent, false, newparent->upper(), equal_to_p);
-
-      CGAL_assertion(false); // Point was not found
-    }
-
-    // Actual removal
-    Leaf_node_handle lnode = static_cast<Leaf_node_handle>(node);
-    if (lnode->size() > 1) {
-      iterator pi = std::find_if(lnode->begin(), lnode->end(), equal_to_p);
-      // FIXME: we should ensure this never happens
-      if (pi == lnode->end()) return false;
-      iterator lasti = lnode->end() - 1;
-      if (pi != lasti) {
-        // Hack to get a non-const iterator
-        std::iter_swap(pts.begin()+(pi-pts.begin()), pts.begin()+(lasti-pts.begin()));
-      }
-      lnode->drop_last_point();
-    } else if (!equal_to_p(*lnode->begin())) {
-      // FIXME: we should ensure this never happens
-      return false;
-    } else if (grandparent) {
-      Node_handle brother = islower ? parent->upper() : parent->lower();
-      if (parent_islower)
-        grandparent->set_lower(brother);
-      else
-        grandparent->set_upper(brother);
-    } else if (parent) {
-      tree_root = islower ? parent->upper() : parent->lower();
-    } else {
-      clear();
-    }
-    return true;
-  }
-
-public:
-  //For efficiency; reserve the size of the points vectors in advance (if the number of points is already known).
-  void reserve(size_t size)
-  {
-    pts.reserve(size);
-  }
-
-  //Get the capacity of the underlying points vector.
-  size_t capacity()
-  {
-    return pts.capacity();
-  }
-
-
-  template <class OutputIterator, class FuzzyQueryItem>
-  OutputIterator
-  search(OutputIterator it, const FuzzyQueryItem& q) const
-  {
-    if(! pts.empty()){
-
-      if(! is_built()){
-        const_build();
-      }
-      Kd_tree_rectangle<FT,D> b(*bbox);
-      return tree_root->search(it,q,b);
-    }
-    return it;
-  }
-
-
-  template <class FuzzyQueryItem>
-  boost::optional<Point_d>
-  search_any_point(const FuzzyQueryItem& q) const
-  {
-    if(! pts.empty()){
-
-      if(! is_built()){
-        const_build();
-      }
-      Kd_tree_rectangle<FT,D> b(*bbox);
-      return tree_root->search_any_point(q,b);
-    }
-    return boost::none;
-  }
-
-
-  ~Kd_tree() {
-    if(is_built()){
-      delete bbox;
-    }
-  }
-
-
-  const SearchTraits&
-  traits() const
-  {
-    return traits_;
-  }
-
-  Node_const_handle
-  root() const
-  {
-    if(! is_built()){
-      const_build();
-    }
-    return tree_root;
-  }
-
-  Node_handle
-  root()
-  {
-    if(! is_built()){
-      build();
-    }
-    return tree_root;
-  }
-
-  void
-  print() const
-  {
-    if(! is_built()){
-      const_build();
-    }
-    root()->print();
-  }
-
-  const Kd_tree_rectangle<FT,D>&
-  bounding_box() const
-  {
-    if(! is_built()){
-      const_build();
-    }
-    return *bbox;
-  }
-
-  const_iterator
-  begin() const
-  {
-    return pts.begin();
-  }
-
-  const_iterator
-  end() const
-  {
-    return pts.end();
-  }
-
-  size_type
-  size() const
-  {
-    return pts.size();
-  }
-
-  // Print statistics of the tree.
-  std::ostream&
-  statistics(std::ostream& s) const
-  {
-    if(! is_built()){
-      const_build();
-    }
-    s << "Tree statistics:" << std::endl;
-    s << "Number of items stored: "
-      << root()->num_items() << std::endl;
-    s << "Number of nodes: "
-      << root()->num_nodes() << std::endl;
-    s << " Tree depth: " << root()->depth() << std::endl;
-    return s;
-  }
-
-
-};
-
-} // namespace CGAL
-
-#endif // CGAL_KD_TREE_H
diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h
deleted file mode 100644
index 909ee260..00000000
--- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h
+++ /dev/null
@@ -1,586 +0,0 @@
-// Copyright (c) 2002,2011  Utrecht University (The Netherlands).
-// All rights reserved.
-//
-// This file is part of CGAL (www.cgal.org).
-// 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.
-//
-// Licensees holding a valid commercial license may use this file in
-// accordance with the commercial license agreement provided with the software.
-//
-// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
-// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
-//
-// $URL$
-// $Id$
-// 
-//
-// Authors       : Hans Tangelder (<hanst@cs.uu.nl>)
-
-#ifndef CGAL_KD_TREE_NODE_H
-#define CGAL_KD_TREE_NODE_H
-
-#include "CGAL/Splitters.h"
-
-#include <CGAL/Compact_container.h>
-#include <boost/cstdint.hpp>
-
-namespace CGAL {
-
-  template <class SearchTraits, class Splitter, class UseExtendedNode>
-  class Kd_tree;
-
-  template < class TreeTraits, class Splitter, class UseExtendedNode >
-  class Kd_tree_node {
-
-     friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>;
-
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_handle Node_handle;
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_const_handle Node_const_handle;
-     typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Internal_node_handle Internal_node_handle;
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Internal_node_const_handle Internal_node_const_handle;
-     typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Leaf_node_handle Leaf_node_handle;
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Leaf_node_const_handle Leaf_node_const_handle;
-    typedef typename TreeTraits::Point_d Point_d;
-
-    typedef typename TreeTraits::FT FT;
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator;
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Point_d_iterator Point_d_iterator;
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::iterator iterator;
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::D D;
-
-    bool leaf;
-
-  public :
-    Kd_tree_node(bool leaf_)
-      :leaf(leaf_){}
-
-    bool is_leaf() const{
-      return leaf;
-    }
-
-    std::size_t
-    num_items() const
-    {
-      if (is_leaf()){
-        Leaf_node_const_handle node =
-          static_cast<Leaf_node_const_handle>(this);
-        return node->size();
-      }
-      else {
-        Internal_node_const_handle node =
-          static_cast<Internal_node_const_handle>(this);
-        return node->lower()->num_items() + node->upper()->num_items();
-      }
-    }
-
-    std::size_t
-    num_nodes() const
-    {
-      if (is_leaf()) return 1;
-      else {
-        Internal_node_const_handle node =
-          static_cast<Internal_node_const_handle>(this);
-        return node->lower()->num_nodes() + node->upper()->num_nodes();
-      }
-    }
-
-    int
-    depth(const int current_max_depth) const
-    {
-      if (is_leaf()){
-        return current_max_depth;
-      }
-      else {
-        Internal_node_const_handle node =
-          static_cast<Internal_node_const_handle>(this);
-        return
-             (std::max)( node->lower()->depth(current_max_depth + 1),
-                         node->upper()->depth(current_max_depth + 1));
-      }
-    }
-
-    int
-    depth() const
-    {
-      return depth(1);
-    }
-
-    template <class OutputIterator>
-    OutputIterator
-    tree_items(OutputIterator it) const {
-      if (is_leaf()) {
-         Leaf_node_const_handle node =
-          static_cast<Leaf_node_const_handle>(this);
-         if (node->size()>0)
-            for (iterator i=node->begin(); i != node->end(); i++)
-              {*it=*i; ++it;}
-        }
-      else {
-         Internal_node_const_handle node =
-          static_cast<Internal_node_const_handle>(this);
-          it=node->lower()->tree_items(it);
-          it=node->upper()->tree_items(it);
-      }
-      return it;
-    }
-
-
-    boost::optional<Point_d>
-    any_tree_item() const {
-      boost::optional<Point_d> result = boost::none;
-      if (is_leaf()) {
-         Leaf_node_const_handle node =
-          static_cast<Leaf_node_const_handle>(this);
-         if (node->size()>0){
-           return boost::make_optional(*(node->begin()));
-         }
-        }
-      else {
-         Internal_node_const_handle node =
-          static_cast<Internal_node_const_handle>(this);
-          result = node->lower()->any_tree_item();
-          if(! result){
-            result = node->upper()->any_tree_item();
-          }
-      }
-      return result;
-    }
-
-
-     void
-    indent(int d) const
-    {
-      for(int i = 0; i < d; i++){
-        std::cout << " ";
-      }
-    }
-
-
-    void
-    print(int d = 0) const
-    {
-      if (is_leaf()) {
-        Leaf_node_const_handle node =
-          static_cast<Leaf_node_const_handle>(this);
-        indent(d);
-        std::cout << "leaf" << std::endl;
-        if (node->size()>0)
-          for (iterator i=node->begin(); i != node->end(); i++)
-          {indent(d);std::cout << *i << std::endl;}
-      }
-      else {
-        Internal_node_const_handle node =
-          static_cast<Internal_node_const_handle>(this);
-        indent(d);
-        std::cout << "lower tree" << std::endl;
-        node->lower()->print(d+1);
-        indent(d);
-        std::cout << "separator: dim = " << node->cutting_dimension() << "  val = " << node->cutting_value() << std::endl;
-        indent(d);
-        std::cout << "upper tree" << std::endl;
-        node->upper()->print(d+1);
-      }
-    }
-
-
-    template <class OutputIterator, class FuzzyQueryItem>
-    OutputIterator
-    search(OutputIterator it, const FuzzyQueryItem& q,
-           Kd_tree_rectangle<FT,D>& b) const
-    {
-      if (is_leaf()) {
-        Leaf_node_const_handle node =
-          static_cast<Leaf_node_const_handle>(this);
-        if (node->size()>0)
-          for (iterator i=node->begin(); i != node->end(); i++)
-            if (q.contains(*i))
-              {*it++=*i;}
-      }
-      else {
-         Internal_node_const_handle node =
-          static_cast<Internal_node_const_handle>(this);
-        // after splitting b denotes the lower part of b
-        Kd_tree_rectangle<FT,D> b_upper(b);
-        b.split(b_upper, node->cutting_dimension(),
-                node->cutting_value());
-
-        if (q.outer_range_contains(b))
-          it=node->lower()->tree_items(it);
-        else
-          if (q.inner_range_intersects(b))
-            it=node->lower()->search(it,q,b);
-        if  (q.outer_range_contains(b_upper))
-          it=node->upper()->tree_items(it);
-        else
-          if (q.inner_range_intersects(b_upper))
-            it=node->upper()->search(it,q,b_upper);
-      };
-      return it;
-    }
-
-
-    template <class FuzzyQueryItem>
-    boost::optional<Point_d>
-    search_any_point(const FuzzyQueryItem& q,
-                     Kd_tree_rectangle<FT,D>& b) const
-    {
-      boost::optional<Point_d> result = boost::none;
-      if (is_leaf()) {
-        Leaf_node_const_handle node =
-          static_cast<Leaf_node_const_handle>(this);
-        if (node->size()>0)
-          for (iterator i=node->begin(); i != node->end(); i++)
-            if (q.contains(*i))
-              { result = *i; break; }
-      }
-      else {
-         Internal_node_const_handle node =
-          static_cast<Internal_node_const_handle>(this);
-        // after splitting b denotes the lower part of b
-        Kd_tree_rectangle<FT,D> b_upper(b);
-        b.split(b_upper, node->cutting_dimension(),
-                node->cutting_value());
-
-        if (q.outer_range_contains(b)){
-          result = node->lower()->any_tree_item();
-        }else{
-          if (q.inner_range_intersects(b)){
-            result = node->lower()->search_any_point(q,b);
-          }
-        }
-        if(result){
-          return result;
-        }
-        if  (q.outer_range_contains(b_upper)){
-          result = node->upper()->any_tree_item();
-        }else{
-          if (q.inner_range_intersects(b_upper))
-            result = node->upper()->search_any_point(q,b_upper);
-        }
-      }
-      return result;
-    }
-
-  };
-
-
-  template < class TreeTraits, class Splitter, class UseExtendedNode >
-  class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{
-
-    friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>;
-
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::iterator iterator;
-    typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base;
-    typedef typename TreeTraits::Point_d Point_d;
-
-  private:
-
-    // private variables for leaf nodes
-    boost::int32_t n; // denotes number of items in a leaf node
-    iterator data; // iterator to data in leaf node
-
-
-  public:
-
-    // default constructor
-    Kd_tree_leaf_node()
-    {}
-
-    Kd_tree_leaf_node(bool leaf_ )
-      : Base(leaf_)
-    {}
-
-    Kd_tree_leaf_node(bool leaf_,unsigned int n_ )
-      : Base(leaf_), n(n_)
-    {}
-
-    // members for all nodes
-
-    // members for leaf nodes only
-    inline
-    unsigned int
-    size() const
-    {
-      return n;
-    }
-
-    inline
-    iterator
-    begin() const
-    {
-      return data;
-    }
-
-    inline
-    iterator
-    end() const
-    {
-      return data + n;
-    }
-
-    inline
-    void
-    drop_last_point()
-    {
-      --n;
-    }
-
-  }; //leaf node
-
-
-
-  template < class TreeTraits, class Splitter, class UseExtendedNode>
-  class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{
-
-    friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>;
-
-    typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base;
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_handle Node_handle;
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_const_handle Node_const_handle;
-
-    typedef typename TreeTraits::FT FT;
-    typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator;
-
-  private:
-
-       // private variables for internal nodes
-    boost::int32_t cut_dim;
-    FT cut_val;
-    Node_handle lower_ch, upper_ch;
-
-
-    // private variables for extended internal nodes
-    FT upper_low_val;
-    FT upper_high_val;
-    FT lower_low_val;
-    FT lower_high_val;
-
-
-  public:
-
-    // default constructor
-    Kd_tree_internal_node()
-    {}
-
-    Kd_tree_internal_node(bool leaf_)
-      : Base(leaf_)
-    {}
-
-
-    // members for internal node and extended internal node
-
-    inline
-    Node_const_handle
-    lower() const
-    {
-      return lower_ch;
-    }
-
-    inline
-    Node_const_handle
-    upper() const
-    {
-      return upper_ch;
-    }
-
-    inline
-    Node_handle
-    lower()
-    {
-      return lower_ch;
-    }
-
-    inline
-    Node_handle
-    upper()
-    {
-      return upper_ch;
-    }
-
-    inline
-    void
-    set_lower(Node_handle nh)
-    {
-      lower_ch = nh;
-    }
-
-    inline
-    void
-    set_upper(Node_handle nh)
-    {
-      upper_ch = nh;
-    }
-
-    // inline Separator& separator() {return sep; }
-    // use instead
-    inline
-    void set_separator(Separator& sep){
-      cut_dim = sep.cutting_dimension();
-      cut_val = sep.cutting_value();
-    }
-
-    inline
-    FT
-    cutting_value() const
-    {
-      return cut_val;
-    }
-
-    inline
-    int
-    cutting_dimension() const
-    {
-      return cut_dim;
-    }
-
-    // members for extended internal node only
-    inline
-    FT
-    upper_low_value() const
-    {
-      return upper_low_val;
-    }
-
-    inline
-    FT
-    upper_high_value() const
-    {
-      return upper_high_val;
-    }
-
-    inline
-    FT
-    lower_low_value() const
-    {
-      return lower_low_val;
-    }
-
-    inline
-    FT
-    lower_high_value() const
-    {
-      return lower_high_val;
-    }
-
-    /*Separator&
-    separator()
-    {
-      return Separator(cutting_dimension,cutting_value);
-    }*/
-
-
-  };//internal node
-
- template < class TreeTraits, class Splitter>
- class Kd_tree_internal_node<TreeTraits,Splitter,Tag_false> : public Kd_tree_node< TreeTraits, Splitter, Tag_false >{
-
-    friend class Kd_tree<TreeTraits,Splitter,Tag_false>;
-
-    typedef Kd_tree_node< TreeTraits, Splitter, Tag_false> Base;
-    typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Node_handle Node_handle;
-    typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Node_const_handle Node_const_handle;
-
-    typedef typename TreeTraits::FT FT;
-    typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Separator Separator;
-
-  private:
-
-       // private variables for internal nodes
-    boost::uint8_t cut_dim;
-    FT cut_val;
-
-    Node_handle lower_ch, upper_ch;
-
-  public:
-
-    // default constructor
-    Kd_tree_internal_node()
-    {}
-
-    Kd_tree_internal_node(bool leaf_)
-      : Base(leaf_)
-    {}
-
-
-    // members for internal node and extended internal node
-
-    inline
-    Node_const_handle
-    lower() const
-    {
-      return lower_ch;
-    }
-
-    inline
-    Node_const_handle
-    upper() const
-    {
-      return upper_ch;
-    }
-
-    inline
-    Node_handle
-    lower()
-    {
-      return lower_ch;
-    }
-
-    inline
-    Node_handle
-    upper()
-    {
-      return upper_ch;
-    }
-
-    inline
-    void
-    set_lower(Node_handle nh)
-    {
-      lower_ch = nh;
-    }
-
-    inline
-    void
-    set_upper(Node_handle nh)
-    {
-      upper_ch = nh;
-    }
-
-    // inline Separator& separator() {return sep; }
-    // use instead
-
-    inline
-    void set_separator(Separator& sep){
-      cut_dim = sep.cutting_dimension();
-      cut_val = sep.cutting_value();
-    }
-
-    inline
-    FT
-    cutting_value() const
-    {
-      return cut_val;
-    }
-
-    inline
-    int
-    cutting_dimension() const
-    {
-      return cut_dim;
-    }
-
-   /* Separator&
-    separator()
-    {
-      return Separator(cutting_dimension,cutting_value);
-    }*/
-
-
-  };//internal node
-
-
-
-} // namespace CGAL
-#endif // CGAL_KDTREE_NODE_H
diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h b/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h
deleted file mode 100644
index dbe707ed..00000000
--- a/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h
+++ /dev/null
@@ -1,621 +0,0 @@
-// Copyright (c) 2002,2011 Utrecht University (The Netherlands).
-// All rights reserved.
-//
-// This file is part of CGAL (www.cgal.org).
-// 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.
-//
-// Licensees holding a valid commercial license may use this file in
-// accordance with the commercial license agreement provided with the software.
-//
-// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
-// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
-//
-// $URL$
-// $Id$
-// 
-//
-// Author(s)     : Hans Tangelder (<hanst@cs.uu.nl>)
-
-#ifndef CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH
-#define CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH
-
-#include "CGAL/Kd_tree.h"
-
-#include <cstring>
-#include <list>
-#include <queue>
-#include <memory>
-#include <CGAL/Euclidean_distance.h>
-#include <CGAL/tuple.h>
-
-namespace CGAL {
-
-  template <class SearchTraits,
-            class Distance_= typename internal::Spatial_searching_default_distance<SearchTraits>::type,
-            class Splitter_ = Sliding_midpoint<SearchTraits>,
-            class Tree_= Kd_tree<SearchTraits, Splitter_, Tag_true> >
-  class Orthogonal_incremental_neighbor_search {
-
-  public:
-    typedef Splitter_ Splitter;
-    typedef Tree_  Tree;
-    typedef Distance_ Distance;
-    typedef typename SearchTraits::Point_d Point_d;
-    typedef typename Distance::Query_item Query_item;
-    typedef typename SearchTraits::FT FT;
-    typedef typename Tree::Point_d_iterator Point_d_iterator;
-    typedef typename Tree::Node_const_handle Node_const_handle;
-
-    typedef std::pair<Point_d,FT> Point_with_transformed_distance;
-    typedef CGAL::cpp11::tuple<Node_const_handle,FT,std::vector<FT> > Node_with_distance;
-    typedef std::vector<Node_with_distance*> Node_with_distance_vector;
-    typedef std::vector<Point_with_transformed_distance*> Point_with_transformed_distance_vector;
-
-    template<class T>
-    struct Object_wrapper
-    {
-      T object;
-      Object_wrapper(const T& t):object(t){}
-      const T& operator* () const { return object; }
-      const T* operator-> () const { return &object; }
-    };
-
-    class Iterator_implementation {
-      SearchTraits traits;
-    public:
-
-      int number_of_neighbours_computed;
-      int number_of_internal_nodes_visited;
-      int number_of_leaf_nodes_visited;
-      int number_of_items_visited;
-
-    private:
-
-      typedef std::vector<FT> Distance_vector;
-
-      Distance_vector dists;
-
-      Distance Orthogonal_distance_instance;
-
-      FT multiplication_factor;
-
-      Query_item query_point;
-
-      FT distance_to_root;
-
-      bool search_nearest_neighbour;
-
-      FT rd;
-
-
-      class Priority_higher {
-      public:
-
-        bool search_nearest;
-
-        Priority_higher(bool search_the_nearest_neighbour)
-          : search_nearest(search_the_nearest_neighbour)
-        {}
-
-        //highest priority is smallest distance
-        bool
-        operator() (Node_with_distance* n1, Node_with_distance* n2) const
-        {
-          return (search_nearest) ? (CGAL::cpp11::get<1>(*n1) > CGAL::cpp11::get<1>(*n2)) : (CGAL::cpp11::get<1>(*n2) > CGAL::cpp11::get<1>(*n1));
-        }
-      };
-
-      class Distance_smaller {
-
-      public:
-
-        bool search_nearest;
-
-        Distance_smaller(bool search_the_nearest_neighbour)
-          : search_nearest(search_the_nearest_neighbour)
-        {}
-
-        //highest priority is smallest distance
-        bool operator() (Point_with_transformed_distance* p1, Point_with_transformed_distance* p2) const
-        {
-          return (search_nearest) ? (p1->second > p2->second) : (p2->second > p1->second);
-        }
-      };
-
-
-      std::priority_queue<Node_with_distance*, Node_with_distance_vector,
-                          Priority_higher> PriorityQueue;
-
-    public:
-      std::priority_queue<Point_with_transformed_distance*, Point_with_transformed_distance_vector,
-                          Distance_smaller> Item_PriorityQueue;
-
-
-    public:
-
-      int reference_count;
-
-
-
-      // constructor
-      Iterator_implementation(const Tree& tree,const Query_item& q, const Distance& tr,
-                              FT Eps=FT(0.0), bool search_nearest=true)
-        : traits(tree.traits()),number_of_neighbours_computed(0), number_of_internal_nodes_visited(0),
-        number_of_leaf_nodes_visited(0), number_of_items_visited(0),
-        Orthogonal_distance_instance(tr), multiplication_factor(Orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)),
-        query_point(q), search_nearest_neighbour(search_nearest),
-        PriorityQueue(Priority_higher(search_nearest)), Item_PriorityQueue(Distance_smaller(search_nearest)),
-        reference_count(1)
-
-
-      {
-        if (tree.empty()) return;
-
-        typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits.construct_cartesian_const_iterator_d_object();
-        int dim = static_cast<int>(std::distance(ccci(q), ccci(q,0)));
-
-        dists.resize(dim);
-        for(int i=0 ; i<dim ; ++i){
-          dists[i] = 0;
-        }
-
-        if (search_nearest){
-          distance_to_root=
-            Orthogonal_distance_instance.min_distance_to_rectangle(q, tree.bounding_box(),dists);
-          Node_with_distance *The_Root = new Node_with_distance(tree.root(),
-                                                                distance_to_root, dists);
-          PriorityQueue.push(The_Root);
-
-          // rd is the distance of the top of the priority queue to q
-          rd=CGAL::cpp11::get<1>(*The_Root);
-          Compute_the_next_nearest_neighbour();
-        }
-         else{
-           distance_to_root=
-         Orthogonal_distance_instance.max_distance_to_rectangle(q,
-                                                tree.bounding_box(), dists);
-           Node_with_distance *The_Root = new Node_with_distance(tree.root(),
-                                                                 distance_to_root, dists);
-        PriorityQueue.push(The_Root);
-
-        // rd is the distance of the top of the priority queue to q
-        rd=CGAL::cpp11::get<1>(*The_Root);
-        Compute_the_next_furthest_neighbour();
-         }
-
-
-      }
-
-      // * operator
-      const Point_with_transformed_distance&
-      operator* () const
-      {
-        return *(Item_PriorityQueue.top());
-      }
-
-      // prefix operator
-      Iterator_implementation&
-      operator++()
-      {
-        Delete_the_current_item_top();
-        if(search_nearest_neighbour)
-          Compute_the_next_nearest_neighbour();
-        else
-          Compute_the_next_furthest_neighbour();
-        return *this;
-      }
-
-      // postfix operator
-      Object_wrapper<Point_with_transformed_distance>
-      operator++(int)
-      {
-        Object_wrapper<Point_with_transformed_distance> result( *(Item_PriorityQueue.top()) );
-        ++*this;
-        return result;
-      }
-
-      // Print statistics of the general priority search process.
-      std::ostream&
-      statistics (std::ostream& s) const {
-        s << "Orthogonal priority search statistics:"
-          << std::endl;
-        s << "Number of internal nodes visited:"
-          << number_of_internal_nodes_visited << std::endl;
-        s << "Number of leaf nodes visited:"
-          << number_of_leaf_nodes_visited << std::endl;
-        s << "Number of items visited:"
-          << number_of_items_visited << std::endl;
-        s << "Number of neighbours computed:"
-          << number_of_neighbours_computed << std::endl;
-        return s;
-      }
-
-
-      //destructor
-      ~Iterator_implementation()
-      {
-        while (!PriorityQueue.empty()) {
-          Node_with_distance* The_top=PriorityQueue.top();
-          PriorityQueue.pop();
-          delete The_top;
-        }
-        while (!Item_PriorityQueue.empty()) {
-          Point_with_transformed_distance* The_top=Item_PriorityQueue.top();
-          Item_PriorityQueue.pop();
-          delete The_top;
-        }
-      }
-
-    private:
-
-      void
-      Delete_the_current_item_top()
-      {
-        Point_with_transformed_distance* The_item_top=Item_PriorityQueue.top();
-        Item_PriorityQueue.pop();
-        delete The_item_top;
-      }
-
-      void
-      Compute_the_next_nearest_neighbour()
-      {
-        // compute the next item
-        bool next_neighbour_found=false;
-        if (!(Item_PriorityQueue.empty())) {
-            next_neighbour_found=
-              (multiplication_factor*rd > Item_PriorityQueue.top()->second);
-        }
-        typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
-        typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point);
-        // otherwise browse the tree further
-        while ((!next_neighbour_found) && (!PriorityQueue.empty())) {
-          Node_with_distance* The_node_top=PriorityQueue.top();
-          Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top);
-          dists = CGAL::cpp11::get<2>(*The_node_top);
-          PriorityQueue.pop();
-          delete The_node_top;
-          FT copy_rd=rd;
-          while (!(N->is_leaf())) { // compute new distance
-            typename Tree::Internal_node_const_handle node =
-            static_cast<typename Tree::Internal_node_const_handle>(N);
-            number_of_internal_nodes_visited++;
-            int new_cut_dim=node->cutting_dimension();
-            FT new_rd,dst = dists[new_cut_dim];
-            FT val = *(query_point_it + new_cut_dim);
-            FT diff1 = val - node->upper_low_value();
-            FT diff2 = val - node->lower_high_value();
-            if (diff1 + diff2 < FT(0.0)) {
-              new_rd=
-                Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim);
-
-              CGAL_assertion(new_rd >= copy_rd);
-                dists[new_cut_dim] = diff1;
-                Node_with_distance *Upper_Child =
-                  new Node_with_distance(node->upper(), new_rd, dists);
-                PriorityQueue.push(Upper_Child);
-                dists[new_cut_dim] = dst;
-                N=node->lower();
-
-            }
-            else { // compute new distance
-              new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim);
-              CGAL_assertion(new_rd >= copy_rd);
-                dists[new_cut_dim] = diff2;
-                Node_with_distance *Lower_Child =
-                  new Node_with_distance(node->lower(), new_rd, dists);
-                PriorityQueue.push(Lower_Child);
-                dists[new_cut_dim] = dst;
-                N=node->upper();
-            }
-          }
-          // n is a leaf
-          typename Tree::Leaf_node_const_handle node =
-            static_cast<typename Tree::Leaf_node_const_handle>(N);
-          number_of_leaf_nodes_visited++;
-          if (node->size() > 0) {
-            for (typename Tree::iterator it=node->begin(); it != node->end(); it++) {
-              number_of_items_visited++;
-              FT distance_to_query_point=
-                Orthogonal_distance_instance.transformed_distance(query_point,*it);
-              Point_with_transformed_distance *NN_Candidate=
-                new Point_with_transformed_distance(*it,distance_to_query_point);
-              Item_PriorityQueue.push(NN_Candidate);
-            }
-            // old top of PriorityQueue has been processed,
-            // hence update rd
-
-            if (!(PriorityQueue.empty()))  {
-              rd = CGAL::cpp11::get<1>(*PriorityQueue.top());
-                next_neighbour_found =
-                  (multiplication_factor*rd >
-                   Item_PriorityQueue.top()->second);
-            }
-            else // priority queue empty => last neighbour found
-              {
-                next_neighbour_found=true;
-              }
-
-            number_of_neighbours_computed++;
-          }
-        }   // next_neighbour_found or priority queue is empty
-        // in the latter case also the item priority quee is empty
-      }
-
-
-      void
-      Compute_the_next_furthest_neighbour()
-      {
-        // compute the next item
-        bool next_neighbour_found=false;
-        if (!(Item_PriorityQueue.empty())) {
-            next_neighbour_found=
-              (rd < multiplication_factor*Item_PriorityQueue.top()->second);
-        }
-        typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
-        typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point);
-        // otherwise browse the tree further
-        while ((!next_neighbour_found) && (!PriorityQueue.empty())) {
-          Node_with_distance* The_node_top=PriorityQueue.top();
-          Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top);
-          dists = CGAL::cpp11::get<2>(*The_node_top);
-          PriorityQueue.pop();
-          delete The_node_top;
-          FT copy_rd=rd;
-          while (!(N->is_leaf())) { // compute new distance
-            typename Tree::Internal_node_const_handle node =
-              static_cast<typename Tree::Internal_node_const_handle>(N);
-            number_of_internal_nodes_visited++;
-            int new_cut_dim=node->cutting_dimension();
-            FT new_rd,dst = dists[new_cut_dim];
-            FT val = *(query_point_it + new_cut_dim);
-            FT diff1 = val - node->upper_low_value();
-            FT diff2 = val - node->lower_high_value();
-            if (diff1 + diff2 < FT(0.0)) {
-              diff1 = val - node->upper_high_value();
-              new_rd=
-                Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim);
-                Node_with_distance *Lower_Child =
-                  new Node_with_distance(node->lower(), copy_rd, dists);
-                PriorityQueue.push(Lower_Child);
-                N=node->upper();
-                dists[new_cut_dim] = diff1;
-                copy_rd=new_rd;
-
-            }
-            else { // compute new distance
-              diff2 = val - node->lower_low_value();
-              new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim);
-                Node_with_distance *Upper_Child =
-                  new Node_with_distance(node->upper(), copy_rd, dists);
-                PriorityQueue.push(Upper_Child);
-                N=node->lower();
-                dists[new_cut_dim] = diff2;
-                copy_rd=new_rd;
-            }
-          }
-          // n is a leaf
-          typename Tree::Leaf_node_const_handle node =
-            static_cast<typename Tree::Leaf_node_const_handle>(N);
-          number_of_leaf_nodes_visited++;
-          if (node->size() > 0) {
-            for (typename Tree::iterator it=node->begin(); it != node->end(); it++) {
-              number_of_items_visited++;
-              FT distance_to_query_point=
-                Orthogonal_distance_instance.transformed_distance(query_point,*it);
-              Point_with_transformed_distance *NN_Candidate=
-                new Point_with_transformed_distance(*it,distance_to_query_point);
-              Item_PriorityQueue.push(NN_Candidate);
-            }
-            // old top of PriorityQueue has been processed,
-            // hence update rd
-
-            if (!(PriorityQueue.empty()))  {
-              rd = CGAL::cpp11::get<1>(*PriorityQueue.top());
-                next_neighbour_found =
-                  (multiplication_factor*rd <
-                   Item_PriorityQueue.top()->second);
-            }
-            else // priority queue empty => last neighbour found
-              {
-                next_neighbour_found=true;
-              }
-
-            number_of_neighbours_computed++;
-          }
-        }   // next_neighbour_found or priority queue is empty
-        // in the latter case also the item priority quee is empty
-      }
-    }; // class Iterator_implementaion
-
-
-
-
-
-
-
-
-
-  public:
-    class iterator;
-    typedef iterator const_iterator;
-
-    // constructor
-    Orthogonal_incremental_neighbor_search(const Tree& tree,
-                                           const Query_item& q, FT Eps = FT(0.0),
-                                           bool search_nearest=true, const Distance& tr=Distance())
-      : m_tree(tree),m_query(q),m_dist(tr),m_Eps(Eps),m_search_nearest(search_nearest)
-    {}
-
-    iterator
-    begin() const
-    {
-      return iterator(m_tree,m_query,m_dist,m_Eps,m_search_nearest);
-    }
-
-    iterator
-    end() const
-    {
-      return iterator();
-    }
-
-    std::ostream&
-    statistics(std::ostream& s)
-    {
-      begin()->statistics(s);
-      return s;
-    }
-
-
-
-
-    class iterator {
-
-    public:
-
-      typedef std::input_iterator_tag iterator_category;
-      typedef Point_with_transformed_distance       value_type;
-      typedef Point_with_transformed_distance*      pointer;
-      typedef const Point_with_transformed_distance&      reference;
-      typedef std::size_t               size_type;
-      typedef std::ptrdiff_t            difference_type;
-      typedef int distance_type;
-
-      //class Iterator_implementation;
-      Iterator_implementation *Ptr_implementation;
-
-
-    public:
-
-      // default constructor
-      iterator()
-      : Ptr_implementation(0)
-      {}
-
-      int
-      the_number_of_items_visited()
-      {
-        return Ptr_implementation->number_of_items_visited;
-      }
-
-      // constructor
-      iterator(const Tree& tree,const Query_item& q, const Distance& tr=Distance(), FT eps=FT(0.0),
-               bool search_nearest=true)
-        : Ptr_implementation(new Iterator_implementation(tree, q, tr, eps, search_nearest))
-        {}
-
-      // copy constructor
-      iterator(const iterator& Iter)
-      {
-        Ptr_implementation = Iter.Ptr_implementation;
-        if (Ptr_implementation != 0) Ptr_implementation->reference_count++;
-      }
-
-      iterator& operator=(const iterator& Iter)
-      {
-        if (Ptr_implementation != Iter.Ptr_implementation){
-          if (Ptr_implementation != 0 && --(Ptr_implementation->reference_count)==0) {
-              delete Ptr_implementation;
-          }
-          Ptr_implementation = Iter.Ptr_implementation;
-          if (Ptr_implementation != 0) Ptr_implementation->reference_count++;
-        }
-        return *this;
-      }
-
-
-      const Point_with_transformed_distance&
-      operator* () const
-      {
-        return *(*Ptr_implementation);
-      }
-
-      // -> operator
-      const Point_with_transformed_distance*
-      operator-> () const
-      {
-        return &*(*Ptr_implementation);
-      }
-
-      // prefix operator
-      iterator&
-      operator++()
-      {
-        ++(*Ptr_implementation);
-        return *this;
-      }
-
-      // postfix operator
-      Object_wrapper<Point_with_transformed_distance>
-      operator++(int)
-      {
-        return (*Ptr_implementation)++;
-      }
-
-
-      bool
-      operator==(const iterator& It) const
-      {
-        if (
-            ((Ptr_implementation == 0) ||
-             Ptr_implementation->Item_PriorityQueue.empty()) &&
-            ((It.Ptr_implementation == 0) ||
-             It.Ptr_implementation->Item_PriorityQueue.empty())
-            )
-          return true;
-        // else
-        return (Ptr_implementation == It.Ptr_implementation);
-      }
-
-      bool
-      operator!=(const iterator& It) const
-      {
-        return !(*this == It);
-      }
-
-      std::ostream&
-      statistics (std::ostream& s)
-      {
-        Ptr_implementation->statistics(s);
-        return s;
-      }
-
-      ~iterator()
-      {
-        if (Ptr_implementation != 0) {
-          Ptr_implementation->reference_count--;
-          if (Ptr_implementation->reference_count==0) {
-            delete Ptr_implementation;
-            Ptr_implementation = 0;
-          }
-        }
-      }
-
-
-    }; // class iterator
-
-    //data members
-    const Tree& m_tree;
-    Query_item m_query;
-    Distance m_dist;
-    FT m_Eps;
-    bool m_search_nearest;
-  }; // class
-
-  template <class Traits, class Query_item, class Distance>
-  void swap (typename Orthogonal_incremental_neighbor_search<Traits,
-             Query_item, Distance>::iterator& x,
-             typename Orthogonal_incremental_neighbor_search<Traits,
-             Query_item, Distance>::iterator& y)
-  {
-    typename Orthogonal_incremental_neighbor_search<Traits,
-      Query_item, Distance>::iterator::Iterator_implementation
-      *tmp = x.Ptr_implementation;
-    x.Ptr_implementation  = y.Ptr_implementation;
-    y.Ptr_implementation = tmp;
-  }
-
-} // namespace CGAL
-
-#endif // CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH_H
diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h
index e9f455d7..253c89b4 100644
--- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h
+++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h
@@ -4,7 +4,7 @@
  *
  *    Author:       Francois Godi
  *
- *    Copyright (C) 2015  INRIA (France)
+ *    Copyright (C) 2015  INRIA
  *
  *    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
@@ -34,141 +34,141 @@ namespace persistence_diagram {
  * \ingroup bottleneck_distance
  */
 class Graph_matching {
-public:
-    /** \internal \brief Constructor constructing an empty matching. */
-    explicit Graph_matching(Persistence_graph &g);
-    /** \internal \brief Copy operator. */
-    Graph_matching& operator=(const Graph_matching& m);
-    /** \internal \brief Is the matching perfect ? */
-    bool perfect() const;
-    /** \internal \brief Augments the matching with a maximal set of edge-disjoint shortest augmenting paths. */
-    bool multi_augment();
-    /** \internal \brief Sets the maximum length of the edges allowed to be added in the matching, 0 initially. */
-    void set_r(double r);
-
-private:
-    Persistence_graph& g;
-    double r;
-    /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */
-    std::vector<int> v_to_u;
-    /** \internal \brief All the unmatched points in U. */
-    std::list<int> unmatched_in_u;
-
-    /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */
-    Layered_neighbors_finder layering() const;
-    /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */
-    bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth);
-    /** \internal \brief Update the matching with the simple augmenting path given as parameter. */
-    void update(std::vector<int> & path);
+ public:
+  /** \internal \brief Constructor constructing an empty matching. */
+  explicit Graph_matching(Persistence_graph &g);
+  /** \internal \brief Copy operator. */
+  Graph_matching& operator=(const Graph_matching& m);
+  /** \internal \brief Is the matching perfect ? */
+  bool perfect() const;
+  /** \internal \brief Augments the matching with a maximal set of edge-disjoint shortest augmenting paths. */
+  bool multi_augment();
+  /** \internal \brief Sets the maximum length of the edges allowed to be added in the matching, 0 initially. */
+  void set_r(double r);
+
+ private:
+  Persistence_graph& g;
+  double r;
+  /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */
+  std::vector<int> v_to_u;
+  /** \internal \brief All the unmatched points in U. */
+  std::list<int> unmatched_in_u;
+
+  /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */
+  Layered_neighbors_finder layering() const;
+  /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */
+  bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth);
+  /** \internal \brief Update the matching with the simple augmenting path given as parameter. */
+  void update(std::vector<int> & path);
 };
 
 inline Graph_matching::Graph_matching(Persistence_graph& g)
     : g(g), r(0.), v_to_u(g.size(), null_point_index()), 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);
+  for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index)
+    unmatched_in_u.emplace_back(u_point_index);
 }
 
 inline Graph_matching& Graph_matching::operator=(const Graph_matching& m) {
-    g = m.g;
-    r = m.r;
-    v_to_u = m.v_to_u;
-    unmatched_in_u = m.unmatched_in_u;
-    return *this;
+  g = m.g;
+  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();
+  return unmatched_in_u.empty();
 }
 
 inline bool Graph_matching::multi_augment() {
-    if (perfect())
-        return false;
-    Layered_neighbors_finder layered_nf(layering());
-    int max_depth = layered_nf.vlayers_number()*2 - 1;
-    double rn = sqrt(g.size());
-    // verification of a necessary criterion in order to shortcut if possible
-    if (max_depth <0 || (unmatched_in_u.size() > rn && max_depth >= rn))
-        return false;
-    bool successful = false;
-    std::list<int> tries(unmatched_in_u);
-    for (auto it = tries.cbegin(); it != tries.cend(); it++)
-        // 'augment' has side-effects which have to be always executed, don't change order
-        successful =  augment(layered_nf, *it, max_depth) || successful;
-    return successful;
+  if (perfect())
+    return false;
+  Layered_neighbors_finder layered_nf(layering());
+  int max_depth = layered_nf.vlayers_number()*2 - 1;
+  double rn = sqrt(g.size());
+  // verification of a necessary criterion in order to shortcut if possible
+  if (max_depth < 0 || (unmatched_in_u.size() > rn && max_depth >= rn))
+    return false;
+  bool successful = false;
+  std::list<int> tries(unmatched_in_u);
+  for (auto it = tries.cbegin(); it != tries.cend(); it++)
+    // 'augment' has side-effects which have to be always executed, don't change order
+    successful = augment(layered_nf, *it, max_depth) || successful;
+  return successful;
 }
 
 inline void Graph_matching::set_r(double r) {
-    this->r = r;
+  this->r = r;
 }
 
 inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) {
-    //V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices.
-    std::vector<int> path;
-    path.emplace_back(u_start_index);
-    do {
-        if (static_cast<int>(path.size()) > max_depth) {
-            path.pop_back();
-            path.pop_back();
-        }
-        if (path.empty())
-            return false;
-        path.emplace_back(layered_nf.pull_near(path.back(), static_cast<int>(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());
-    //if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex
-    path.pop_back();
-    update(path);
-    return true;
+  // V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices.
+  std::vector<int> path;
+  path.emplace_back(u_start_index);
+  do {
+    if (static_cast<int> (path.size()) > max_depth) {
+      path.pop_back();
+      path.pop_back();
+    }
+    if (path.empty())
+      return false;
+    path.emplace_back(layered_nf.pull_near(path.back(), static_cast<int> (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());
+  // if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex
+  path.pop_back();
+  update(path);
+  return true;
 }
 
 inline Layered_neighbors_finder Graph_matching::layering() const {
-    std::list<int> u_vertices(unmatched_in_u);
-    std::list<int> v_vertices;
-    Neighbors_finder nf(g, r);
-    for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index)
-        nf.add(v_point_index);
-    Layered_neighbors_finder layered_nf(g, r);
-    for(int layer = 0; !u_vertices.empty(); layer++) {
-        // one layer is one step in the BFS
-        for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) {
-            std::vector<int> u_succ(nf.pull_all_near(*it1));
-            for (auto it2 = u_succ.begin(); it2 != u_succ.end(); ++it2) {
-                layered_nf.add(*it2, layer);
-                v_vertices.emplace_back(*it2);
-            }
-        }
-        // When the above for finishes, we have progress of one half-step (from U to V) in the BFS
-        u_vertices.clear();
-        bool end = false;
-        for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++)
-            if (v_to_u.at(*it) == null_point_index())
-                // we stop when a nearest exposed V vertex (from U exposed vertices) has been found
-                end = true;
-            else
-                u_vertices.emplace_back(v_to_u.at(*it));
-        // When the above for finishes, we have progress of one half-step (from V to U) in the BFS
-        if (end)
-            return layered_nf;
-        v_vertices.clear();
+  std::list<int> u_vertices(unmatched_in_u);
+  std::list<int> v_vertices;
+  Neighbors_finder nf(g, r);
+  for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index)
+    nf.add(v_point_index);
+  Layered_neighbors_finder layered_nf(g, r);
+  for (int layer = 0; !u_vertices.empty(); layer++) {
+    // one layer is one step in the BFS
+    for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) {
+      std::vector<int> u_succ(nf.pull_all_near(*it1));
+      for (auto it2 = u_succ.begin(); it2 != u_succ.end(); ++it2) {
+        layered_nf.add(*it2, layer);
+        v_vertices.emplace_back(*it2);
+      }
     }
-    return layered_nf;
+    // When the above for finishes, we have progress of one half-step (from U to V) in the BFS
+    u_vertices.clear();
+    bool end = false;
+    for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++)
+      if (v_to_u.at(*it) == null_point_index())
+        // we stop when a nearest exposed V vertex (from U exposed vertices) has been found
+        end = true;
+      else
+        u_vertices.emplace_back(v_to_u.at(*it));
+    // When the above for finishes, we have progress of one half-step (from V to U) in the BFS
+    if (end)
+      return layered_nf;
+    v_vertices.clear();
+  }
+  return layered_nf;
 }
 
 inline void Graph_matching::update(std::vector<int>& path) {
-    unmatched_in_u.remove(path.front());
-    for (auto it = path.cbegin(); it != path.cend(); ++it) {
-        // Be careful, the iterator is incremented twice each time
-        int tmp = *it;
-        v_to_u[*(++it)] = tmp;
-    }
+  unmatched_in_u.remove(path.front());
+  for (auto it = path.cbegin(); it != path.cend(); ++it) {
+    // Be careful, the iterator is incremented twice each time
+    int tmp = *it;
+    v_to_u[*(++it)] = tmp;
+  }
 }
 
 
diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h
index 70342d0c..0b2d26fe 100644
--- a/src/Bottleneck_distance/include/gudhi/Internal_point.h
+++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h
@@ -4,7 +4,7 @@
  *
  *    Author:       Francois Godi
  *
- *    Copyright (C) 2015  INRIA (France)
+ *    Copyright (C) 2015  INRIA
  *
  *    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
@@ -32,39 +32,60 @@ int null_point_index();
 
 /** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */
 struct Internal_point {
-    double vec[2];
-    int point_index;
-    Internal_point() {}
-    Internal_point(double x, double y, int p_i) { vec[0]=x; vec[1]=y; point_index = p_i; }
-    double x() const { return vec[ 0 ]; }
-    double y() const { return vec[ 1 ]; }
-    double& x() { return vec[ 0 ]; }
-    double& y() { return vec[ 1 ]; }
-    bool operator==(const Internal_point& p) const
-    {
-        return point_index==p.point_index;
-    }
-    bool  operator!=(const Internal_point& p) const { return !(*this == p); }
+  double vec[2];
+  int point_index;
+
+  Internal_point() { }
+
+  Internal_point(double x, double y, int p_i) {
+    vec[0] = x;
+    vec[1] = y;
+    point_index = p_i;
+  }
+
+  double x() const {
+    return vec[ 0 ];
+  }
+
+  double y() const {
+    return vec[ 1 ];
+  }
+
+  double& x() {
+    return vec[ 0 ];
+  }
+
+  double& y() {
+    return vec[ 1 ];
+  }
+
+  bool operator==(const Internal_point& p) const {
+    return point_index == p.point_index;
+  }
+
+  bool operator!=(const Internal_point& p) const {
+    return !(*this == p);
+  }
 };
 
 inline int null_point_index() {
-    return -1;
+  return -1;
 }
 
 struct Construct_coord_iterator {
-    typedef  const double* result_type;
-    const double* operator()(const Internal_point& p) const
-    { return p.vec; }
-    const double* operator()(const Internal_point& p, int)  const
-    { return p.vec+2; }
+  typedef const double* result_type;
+
+  const double* operator()(const Internal_point& p) const {
+    return p.vec;
+  }
+
+  const double* operator()(const Internal_point& p, int) const {
+    return p.vec + 2;
+  }
 };
 
 }  // namespace persistence_diagram
 
 }  // namespace Gudhi
 
-
-
-
-
 #endif  // INTERNAL_POINT_H_
diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
index 792925b7..96ece360 100644
--- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
+++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
@@ -4,7 +4,7 @@
  *
  *    Author:       Francois Godi
  *
- *    Copyright (C) 2015  INRIA (France)
+ *    Copyright (C) 2015  INRIA
  *
  *    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
@@ -24,10 +24,9 @@
 #define NEIGHBORS_FINDER_H_
 
 // Inclusion order is important for CGAL patch
-#include "CGAL/Kd_tree_node.h"
-#include "CGAL/Kd_tree.h"
-#include "CGAL/Orthogonal_k_neighbor_search.h"
-
+#include <CGAL/Kd_tree.h>
+#include <CGAL/Kd_tree_node.h>
+#include <CGAL/Orthogonal_k_neighbor_search.h>
 #include <CGAL/Weighted_Minkowski_distance.h>
 #include <CGAL/Search_traits.h>
 
@@ -43,123 +42,126 @@ namespace persistence_diagram {
 /** \internal \brief data structure used to find any point (including projections) in V near to a query point from U
  * (which can be a projection).
  *
- * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically removed.
+ * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically
+ * removed.
  *
  * \ingroup bottleneck_distance
  */
 class Neighbors_finder {
-
-    typedef CGAL::Dimension_tag<2> D;
-    typedef CGAL::Search_traits<double, Internal_point, const double*, Construct_coord_iterator, D> Traits;
-    typedef CGAL::Weighted_Minkowski_distance<Traits> Distance;
-    typedef CGAL::Orthogonal_k_neighbor_search<Traits, Distance> K_neighbor_search;
-    typedef K_neighbor_search::Tree Kd_tree;
-
-public:
-    /** \internal \brief Constructor taking the near distance definition as parameter. */
-    Neighbors_finder(const Persistence_graph& g, double r);
-    /** \internal \brief A point added will be possibly pulled. */
-    void add(int v_point_index);
-    /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */
-    int pull_near(int u_point_index);
-    /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */
-    std::vector<int> pull_all_near(int u_point_index);
-
-private:
-    const Persistence_graph& g;
-    const double r;
-    Kd_tree kd_t;
-    std::unordered_set<int> projections_f;
+  typedef CGAL::Dimension_tag<2> D;
+  typedef CGAL::Search_traits<double, Internal_point, const double*, Construct_coord_iterator, D> Traits;
+  typedef CGAL::Weighted_Minkowski_distance<Traits> Distance;
+  typedef CGAL::Orthogonal_k_neighbor_search<Traits, Distance> K_neighbor_search;
+  typedef K_neighbor_search::Tree Kd_tree;
+
+ public:
+  /** \internal \brief Constructor taking the near distance definition as parameter. */
+  Neighbors_finder(const Persistence_graph& g, double r);
+  /** \internal \brief A point added will be possibly pulled. */
+  void add(int v_point_index);
+  /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if
+   * there isn't such a point. */
+  int pull_near(int u_point_index);
+  /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */
+  std::vector<int> pull_all_near(int u_point_index);
+
+ private:
+  const Persistence_graph& g;
+  const double r;
+  Kd_tree kd_t;
+  std::unordered_set<int> projections_f;
 };
 
 /** \internal \brief data structure used to find any point (including projections) in V near to a query point from U
  * (which can be a projection) in a layered graph layer given as parmeter.
  *
- * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically removed.
+ * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically
+ * removed.
  *
  * \ingroup bottleneck_distance
  */
 class Layered_neighbors_finder {
-public:
-    /** \internal \brief Constructor taking the near distance definition as parameter. */
-    Layered_neighbors_finder(const Persistence_graph& g, double r);
-    /** \internal \brief A point added will be possibly pulled. */
-    void add(int v_point_index, int vlayer);
-    /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */
-    int pull_near(int u_point_index, int vlayer);
-    /** \internal \brief Returns the number of layers. */
-    int vlayers_number() const;
-
-private:
-    const Persistence_graph& g;
-    const double r;
-    std::vector<std::unique_ptr<Neighbors_finder>> neighbors_finder;
+ public:
+  /** \internal \brief Constructor taking the near distance definition as parameter. */
+  Layered_neighbors_finder(const Persistence_graph& g, double r);
+  /** \internal \brief A point added will be possibly pulled. */
+  void add(int v_point_index, int vlayer);
+  /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if
+   * there isn't such a point. */
+  int pull_near(int u_point_index, int vlayer);
+  /** \internal \brief Returns the number of layers. */
+  int vlayers_number() const;
+
+ private:
+  const Persistence_graph& g;
+  const double r;
+  std::vector<std::unique_ptr<Neighbors_finder>> neighbors_finder;
 };
 
 inline Neighbors_finder::Neighbors_finder(const Persistence_graph& g, double r) :
     g(g), r(r), kd_t(), 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
-        kd_t.insert(g.get_v_point(v_point_index));
+  if (g.on_the_v_diagonal(v_point_index))
+    projections_f.emplace(v_point_index);
+  else
+    kd_t.insert(g.get_v_point(v_point_index));
 }
 
 inline int Neighbors_finder::pull_near(int u_point_index) {
-    int tmp;
-    int c = g.corresponding_point_in_v(u_point_index);
-    if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()){
-        //Any pair of projection is at distance 0
-        tmp = *projections_f.cbegin();
-        projections_f.erase(tmp);
-    }
-    else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)){
-        //Is the query point near to its projection ?
-        tmp = c;
-        projections_f.erase(tmp);
-    }
-    else{
-        //Is the query point near to a V point in the plane ?
-        Internal_point u_point = g.get_u_point(u_point_index);
-        std::array<double, 2> w = { {1., 1.} };
-         K_neighbor_search search(kd_t, u_point, 1, 0., true, Distance(0, 2, w.begin(), w.end()));
-        auto it = search.begin();
-        if(it==search.end() || g.distance(u_point_index, it->first.point_index) > r)
-            return null_point_index();
-        tmp = it->first.point_index;
-        kd_t.remove(g.get_v_point(tmp));
-    }
-    return tmp;
+  int tmp;
+  int c = g.corresponding_point_in_v(u_point_index);
+  if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()) {
+    // Any pair of projection is at distance 0
+    tmp = *projections_f.cbegin();
+    projections_f.erase(tmp);
+  } else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)) {
+    // Is the query point near to its projection ?
+    tmp = c;
+    projections_f.erase(tmp);
+  } else {
+    // Is the query point near to a V point in the plane ?
+    Internal_point u_point = g.get_u_point(u_point_index);
+    std::array<double, 2> w = {
+      {1., 1.}
+    };
+    K_neighbor_search search(kd_t, u_point, 1, 0., true, Distance(0, 2, w.begin(), w.end()));
+    auto it = search.begin();
+    if (it == search.end() || g.distance(u_point_index, it->first.point_index) > r)
+      return null_point_index();
+    tmp = it->first.point_index;
+    kd_t.remove(g.get_v_point(tmp));
+  }
+  return tmp;
 }
 
 inline std::vector<int> Neighbors_finder::pull_all_near(int u_point_index) {
-    std::vector<int> all_pull;
-    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;
+  std::vector<int> all_pull;
+  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;
 }
 
 inline Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_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(std::unique_ptr<Neighbors_finder>(new Neighbors_finder(g, r)));
-    neighbors_finder.at(vlayer)->add(v_point_index);
+  for (int l = neighbors_finder.size(); l <= vlayer; l++)
+    neighbors_finder.emplace_back(std::unique_ptr<Neighbors_finder>(new 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);
+  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 static_cast<int>(neighbors_finder.size());
+  return static_cast<int> (neighbors_finder.size());
 }
 
 }  // namespace persistence_diagram
diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
index 45a4d586..3a4a5fec 100644
--- a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
+++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
@@ -4,7 +4,7 @@
  *
  *    Author:       Francois Godi
  *
- *    Copyright (C) 2015  INRIA (France)
+ *    Copyright (C) 2015  INRIA
  *
  *    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
@@ -31,147 +31,144 @@ namespace Gudhi {
 
 namespace persistence_diagram {
 
-
 /** \internal \brief Structure representing an euclidean bipartite graph containing
  *  the points from the two persistence diagrams (including the projections).
  *
  * \ingroup bottleneck_distance
  */
 class Persistence_graph {
-public:
-    /** \internal \brief Constructor taking 2 Persistence_Diagrams (concept) as parameters. */
-    template<typename Persistence_diagram1, typename Persistence_diagram2>
-    Persistence_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e);
-    /** \internal \brief Is the given point from U the projection of a point in V ? */
-    bool on_the_u_diagonal(int u_point_index) const;
-    /** \internal \brief Is the given point from V the projection of a point in U ? */
-    bool on_the_v_diagonal(int v_point_index) const;
-    /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */
-    int corresponding_point_in_u(int v_point_index) const;
-    /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */
-    int corresponding_point_in_v(int u_point_index) const;
-    /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */
-    double distance(int u_point_index, int v_point_index) const;
-    /** \internal \brief Returns size = |U| = |V|. */
-    int size() const;
-    /** \internal \brief Is there as many infinite points (alive components) in both diagrams ? */
-    double bottleneck_alive() const;
-    /** \internal \brief Returns the O(n^2) sorted distances between the points. */
-    std::vector<double> sorted_distances() const;
-    /** \internal \brief Returns an upper bound for the diameter of the convex hull of all non infinite points */
-    double diameter_bound() const;
-    /** \internal \brief Returns the corresponding internal point */
-    Internal_point get_u_point(int u_point_index) const;
-    /** \internal \brief Returns the corresponding internal point */
-    Internal_point get_v_point(int v_point_index) const;
-
-private:
-    std::vector<Internal_point> u;
-    std::vector<Internal_point> v;
-    double b_alive;
+ public:
+  /** \internal \brief Constructor taking 2 Persistence_Diagrams (concept) as parameters. */
+  template<typename Persistence_diagram1, typename Persistence_diagram2>
+  Persistence_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e);
+  /** \internal \brief Is the given point from U the projection of a point in V ? */
+  bool on_the_u_diagonal(int u_point_index) const;
+  /** \internal \brief Is the given point from V the projection of a point in U ? */
+  bool on_the_v_diagonal(int v_point_index) const;
+  /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */
+  int corresponding_point_in_u(int v_point_index) const;
+  /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */
+  int corresponding_point_in_v(int u_point_index) const;
+  /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */
+  double distance(int u_point_index, int v_point_index) const;
+  /** \internal \brief Returns size = |U| = |V|. */
+  int size() const;
+  /** \internal \brief Is there as many infinite points (alive components) in both diagrams ? */
+  double bottleneck_alive() const;
+  /** \internal \brief Returns the O(n^2) sorted distances between the points. */
+  std::vector<double> sorted_distances() const;
+  /** \internal \brief Returns an upper bound for the diameter of the convex hull of all non infinite points */
+  double diameter_bound() const;
+  /** \internal \brief Returns the corresponding internal point */
+  Internal_point get_u_point(int u_point_index) const;
+  /** \internal \brief Returns the corresponding internal point */
+  Internal_point get_v_point(int v_point_index) const;
+
+ private:
+  std::vector<Internal_point> u;
+  std::vector<Internal_point> v;
+  double b_alive;
 };
 
 template<typename Persistence_diagram1, typename Persistence_diagram2>
 Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1,
                                      const Persistence_diagram2 &diag2, double e)
-    : u(), v(), b_alive(0.)
-{
-    std::vector<double> u_alive;
-    std::vector<double> v_alive;
-    for (auto it = std::begin(diag1); it != std::end(diag1); ++it){
-        if(std::get<1>(*it) == std::numeric_limits<double>::infinity())
-            u_alive.push_back(std::get<0>(*it));
-        else if (std::get<1>(*it) - std::get<0>(*it) > e)
-            u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size()));
-    }
-    for (auto it = std::begin(diag2); it != std::end(diag2); ++it){
-        if(std::get<1>(*it) == std::numeric_limits<double>::infinity())
-            v_alive.push_back(std::get<0>(*it));
-        else if (std::get<1>(*it) - std::get<0>(*it) > e)
-            v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size()));
-    }
-    if (u.size() < v.size())
-        swap(u, v);
-    std::sort(u_alive.begin(), u_alive.end());
-    std::sort(v_alive.begin(), v_alive.end());
-    if(u_alive.size() != v_alive.size())
-        b_alive = std::numeric_limits<double>::infinity();
-    else for(auto it_u=u_alive.cbegin(), it_v=v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v)
-        b_alive = std::max(b_alive, std::fabs(*it_u - *it_v));
+    : u(), v(), b_alive(0.) {
+  std::vector<double> u_alive;
+  std::vector<double> v_alive;
+  for (auto it = std::begin(diag1); it != std::end(diag1); ++it) {
+    if (std::get<1>(*it) == std::numeric_limits<double>::infinity())
+      u_alive.push_back(std::get<0>(*it));
+    else if (std::get<1>(*it) - std::get<0>(*it) > e)
+      u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size()));
+  }
+  for (auto it = std::begin(diag2); it != std::end(diag2); ++it) {
+    if (std::get<1>(*it) == std::numeric_limits<double>::infinity())
+      v_alive.push_back(std::get<0>(*it));
+    else if (std::get<1>(*it) - std::get<0>(*it) > e)
+      v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size()));
+  }
+  if (u.size() < v.size())
+    swap(u, v);
+  std::sort(u_alive.begin(), u_alive.end());
+  std::sort(v_alive.begin(), v_alive.end());
+  if (u_alive.size() != v_alive.size())
+    b_alive = std::numeric_limits<double>::infinity();
+  else for (auto it_u = u_alive.cbegin(), it_v = v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v)
+      b_alive = std::max(b_alive, std::fabs(*it_u - *it_v));
 }
 
 inline bool Persistence_graph::on_the_u_diagonal(int u_point_index) const {
-    return u_point_index >= static_cast<int> (u.size());
+  return u_point_index >= static_cast<int> (u.size());
 }
 
 inline bool Persistence_graph::on_the_v_diagonal(int v_point_index) const {
-    return v_point_index >= static_cast<int> (v.size());
+  return v_point_index >= static_cast<int> (v.size());
 }
 
 inline int Persistence_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());
+  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_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());
+  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_graph::distance(int u_point_index, int v_point_index) const {
-    if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index))
-        return 0.;
-    Internal_point p_u = get_u_point(u_point_index);
-    Internal_point p_v = get_v_point(v_point_index);
-    return std::max(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y()));
+  if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index))
+    return 0.;
+  Internal_point p_u = get_u_point(u_point_index);
+  Internal_point p_v = get_v_point(v_point_index);
+  return std::max(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y()));
 }
 
 inline int Persistence_graph::size() const {
-    return static_cast<int> (u.size() + v.size());
+  return static_cast<int> (u.size() + v.size());
 }
 
-inline double Persistence_graph::bottleneck_alive() const{
-    return b_alive;
+inline double Persistence_graph::bottleneck_alive() const {
+  return b_alive;
 }
 
 inline std::vector<double> Persistence_graph::sorted_distances() const {
-    std::vector<double> distances;
-    distances.push_back(0.); //for empty diagrams
-    for (int u_point_index = 0; u_point_index < size(); ++u_point_index){
-        distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index)));
-        for (int v_point_index = 0; v_point_index < size(); ++v_point_index)
-            distances.push_back(distance(u_point_index, v_point_index));
-    }
-    std::sort(distances.begin(), distances.end());
-    return distances;
+  std::vector<double> distances;
+  distances.push_back(0.);  // for empty diagrams
+  for (int u_point_index = 0; u_point_index < size(); ++u_point_index) {
+    distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index)));
+    for (int v_point_index = 0; v_point_index < size(); ++v_point_index)
+      distances.push_back(distance(u_point_index, v_point_index));
+  }
+  std::sort(distances.begin(), distances.end());
+  return distances;
 }
 
 inline Internal_point Persistence_graph::get_u_point(int u_point_index) const {
-    if (!on_the_u_diagonal(u_point_index))
-        return u.at(u_point_index);
-    Internal_point projector = v.at(corresponding_point_in_v(u_point_index));
-    double m = (projector.x() + projector.y()) / 2.;
-    return Internal_point(m,m,u_point_index);
+  if (!on_the_u_diagonal(u_point_index))
+    return u.at(u_point_index);
+  Internal_point projector = v.at(corresponding_point_in_v(u_point_index));
+  double m = (projector.x() + projector.y()) / 2.;
+  return Internal_point(m, m, u_point_index);
 }
 
 inline Internal_point Persistence_graph::get_v_point(int v_point_index) const {
-    if (!on_the_v_diagonal(v_point_index))
-        return v.at(v_point_index);
-    Internal_point projector = u.at(corresponding_point_in_u(v_point_index));
-    double m = (projector.x() + projector.y()) / 2.;
-    return Internal_point(m,m,v_point_index);
+  if (!on_the_v_diagonal(v_point_index))
+    return v.at(v_point_index);
+  Internal_point projector = u.at(corresponding_point_in_u(v_point_index));
+  double m = (projector.x() + projector.y()) / 2.;
+  return Internal_point(m, m, v_point_index);
 }
 
 inline double Persistence_graph::diameter_bound() const {
-    double max = 0.;
-    for(auto it = u.cbegin(); it != u.cend(); it++)
-        max = std::max(max, it->y());
-    for(auto it = v.cbegin(); it != v.cend(); it++)
-        max = std::max(max, it->y());
-    return max;
+  double max = 0.;
+  for (auto it = u.cbegin(); it != u.cend(); it++)
+    max = std::max(max, it->y());
+  for (auto it = v.cbegin(); it != v.cend(); it++)
+    max = std::max(max, it->y());
+  return max;
 }
 
-
 }  // namespace persistence_diagram
 
 }  // namespace Gudhi
diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt
index 13213075..a6979d3c 100644
--- a/src/Bottleneck_distance/test/CMakeLists.txt
+++ b/src/Bottleneck_distance/test/CMakeLists.txt
@@ -17,7 +17,6 @@ if(CGAL_FOUND)
   if (NOT CGAL_VERSION VERSION_LESS 4.8.0)
     if (EIGEN3_FOUND)
       add_executable ( bottleneckUT bottleneck_unit_test.cpp )
-      add_executable ( bottleneck_chrono bottleneck_chrono.cpp )
       target_link_libraries(bottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
 
       # Unitary tests
diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
index fba1d369..e39613b3 100644
--- a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
+++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
@@ -4,7 +4,7 @@
  *
  *    Author:       Francois Godi
  *
- *    Copyright (C) 2015  INRIA (France)
+ *    Copyright (C) 2015  INRIA
  *
  *    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
@@ -30,138 +30,138 @@
 
 using namespace Gudhi::persistence_diagram;
 
-int n1 = 81; // a natural number >0
-int n2 = 180; // a natural number >0
-double upper_bound = 406.43; // any real >0
+int n1 = 81;  // a natural number >0
+int n2 = 180;  // a natural number >0
+double upper_bound = 406.43;  // any real >0
 
 
-std::uniform_real_distribution<double> unif(0.,upper_bound);
+std::uniform_real_distribution<double> unif(0., upper_bound);
 std::default_random_engine re;
 std::vector< std::pair<double, double> > v1, v2;
 
-BOOST_AUTO_TEST_CASE(persistence_graph){
-    // Random construction
-    for (int i = 0; i < n1; i++) {
-        double a = unif(re);
-        double b = unif(re);
-        v1.emplace_back(std::min(a,b), std::max(a,b));
-    }
-    for (int i = 0; i < n2; i++) {
-        double a = unif(re);
-        double b = unif(re);
-        v2.emplace_back(std::min(a,b), std::max(a,b));
-    }
-    Persistence_graph g(v1, v2, 0.);
-    std::vector<double> d(g.sorted_distances());
-    //
-    BOOST_CHECK(!g.on_the_u_diagonal(n1-1));
-    BOOST_CHECK(!g.on_the_u_diagonal(n1));
-    BOOST_CHECK(!g.on_the_u_diagonal(n2-1));
-    BOOST_CHECK(g.on_the_u_diagonal(n2));
-    BOOST_CHECK(!g.on_the_v_diagonal(n1-1));
-    BOOST_CHECK(g.on_the_v_diagonal(n1));
-    BOOST_CHECK(g.on_the_v_diagonal(n2-1));
-    BOOST_CHECK(g.on_the_v_diagonal(n2));
-    //
-    BOOST_CHECK(g.corresponding_point_in_u(0)==n2);
-    BOOST_CHECK(g.corresponding_point_in_u(n1)==0);
-    BOOST_CHECK(g.corresponding_point_in_v(0)==n1);
-    BOOST_CHECK(g.corresponding_point_in_v(n2)==0);
-    //
-    BOOST_CHECK(g.size()==(n1+n2));
-    //
-    BOOST_CHECK((int) d.size() == (n1+n2)*(n1+n2) + n1 + n2 + 1);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,0))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n1-1))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n1))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n2-1))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n2))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,(n1+n2)-1))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,0))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n1-1))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n1))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n2-1))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n2))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,(n1+n2)-1))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,0))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n1-1))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n1))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n2-1))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n2))>0);
-    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,(n1+n2)-1))>0);
+BOOST_AUTO_TEST_CASE(persistence_graph) {
+  // Random construction
+  for (int i = 0; i < n1; i++) {
+    double a = unif(re);
+    double b = unif(re);
+    v1.emplace_back(std::min(a, b), std::max(a, b));
+  }
+  for (int i = 0; i < n2; i++) {
+    double a = unif(re);
+    double b = unif(re);
+    v2.emplace_back(std::min(a, b), std::max(a, b));
+  }
+  Persistence_graph g(v1, v2, 0.);
+  std::vector<double> d(g.sorted_distances());
+  //
+  BOOST_CHECK(!g.on_the_u_diagonal(n1 - 1));
+  BOOST_CHECK(!g.on_the_u_diagonal(n1));
+  BOOST_CHECK(!g.on_the_u_diagonal(n2 - 1));
+  BOOST_CHECK(g.on_the_u_diagonal(n2));
+  BOOST_CHECK(!g.on_the_v_diagonal(n1 - 1));
+  BOOST_CHECK(g.on_the_v_diagonal(n1));
+  BOOST_CHECK(g.on_the_v_diagonal(n2 - 1));
+  BOOST_CHECK(g.on_the_v_diagonal(n2));
+  //
+  BOOST_CHECK(g.corresponding_point_in_u(0) == n2);
+  BOOST_CHECK(g.corresponding_point_in_u(n1) == 0);
+  BOOST_CHECK(g.corresponding_point_in_v(0) == n1);
+  BOOST_CHECK(g.corresponding_point_in_v(n2) == 0);
+  //
+  BOOST_CHECK(g.size() == (n1 + n2));
+  //
+  BOOST_CHECK((int) d.size() == (n1 + n2)*(n1 + n2) + n1 + n2 + 1);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, 0)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n1 - 1)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n1)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n2 - 1)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, n2)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0, (n1 + n2) - 1)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, 0)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n1 - 1)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n1)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n2 - 1)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, n2)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1, (n1 + n2) - 1)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, 0)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n1 - 1)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n1)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n2 - 1)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, n2)) > 0);
+  BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1 + n2) - 1, (n1 + n2) - 1)) > 0);
 }
 
 BOOST_AUTO_TEST_CASE(neighbors_finder) {
-    Persistence_graph g(v1, v2, 0.);
-    Neighbors_finder nf(g, 1.);
-    for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2)
-        nf.add(v_point_index);
-    //
-    int v_point_index_1 = nf.pull_near(n2/2);
-    BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.));
-    std::vector<int> l = nf.pull_all_near(n2/2);
-    bool v = true;
-    for(auto it = l.cbegin(); it != l.cend(); ++it)
-        v = v && (g.distance(n2/2,*it)>1.);
-    BOOST_CHECK(v);
-    int v_point_index_2 = nf.pull_near(n2/2);
-    BOOST_CHECK(v_point_index_2 == -1);
+  Persistence_graph g(v1, v2, 0.);
+  Neighbors_finder nf(g, 1.);
+  for (int v_point_index = 1; v_point_index < ((n2 + n1)*9 / 10); v_point_index += 2)
+    nf.add(v_point_index);
+  //
+  int v_point_index_1 = nf.pull_near(n2 / 2);
+  BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.));
+  std::vector<int> l = nf.pull_all_near(n2 / 2);
+  bool v = true;
+  for (auto it = l.cbegin(); it != l.cend(); ++it)
+    v = v && (g.distance(n2 / 2, *it) > 1.);
+  BOOST_CHECK(v);
+  int v_point_index_2 = nf.pull_near(n2 / 2);
+  BOOST_CHECK(v_point_index_2 == -1);
 }
 
 BOOST_AUTO_TEST_CASE(layered_neighbors_finder) {
-    Persistence_graph g(v1, v2, 0.);
-    Layered_neighbors_finder lnf(g, 1.);
-    for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2)
-        lnf.add(v_point_index, v_point_index % 7);
-    //
-    int v_point_index_1 = lnf.pull_near(n2/2,6);
-    BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.));
-    int v_point_index_2 = lnf.pull_near(n2/2,6);
-    BOOST_CHECK(v_point_index_2 == -1);
-    v_point_index_1 = lnf.pull_near(n2/2,0);
-    BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.));
-    v_point_index_2 = lnf.pull_near(n2/2,0);
-    BOOST_CHECK(v_point_index_2 == -1);
+  Persistence_graph g(v1, v2, 0.);
+  Layered_neighbors_finder lnf(g, 1.);
+  for (int v_point_index = 1; v_point_index < ((n2 + n1)*9 / 10); v_point_index += 2)
+    lnf.add(v_point_index, v_point_index % 7);
+  //
+  int v_point_index_1 = lnf.pull_near(n2 / 2, 6);
+  BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.));
+  int v_point_index_2 = lnf.pull_near(n2 / 2, 6);
+  BOOST_CHECK(v_point_index_2 == -1);
+  v_point_index_1 = lnf.pull_near(n2 / 2, 0);
+  BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.));
+  v_point_index_2 = lnf.pull_near(n2 / 2, 0);
+  BOOST_CHECK(v_point_index_2 == -1);
 }
 
 BOOST_AUTO_TEST_CASE(graph_matching) {
-    Persistence_graph g(v1, v2, 0.);
-    Graph_matching m1(g);
-    m1.set_r(0.);
-    int e  = 0;
-    while (m1.multi_augment())
-        ++e;
-    BOOST_CHECK(e > 0);
-    BOOST_CHECK(e <= 2*sqrt(2*(n1+n2)));
-    Graph_matching m2 = m1;
-    BOOST_CHECK(!m2.multi_augment());
-    m2.set_r(upper_bound);
-    e  = 0;
-    while (m2.multi_augment())
-        ++e;
-    BOOST_CHECK(e <= 2*sqrt(2*(n1+n2)));
-    BOOST_CHECK(m2.perfect());
-    BOOST_CHECK(!m1.perfect());
+  Persistence_graph g(v1, v2, 0.);
+  Graph_matching m1(g);
+  m1.set_r(0.);
+  int e = 0;
+  while (m1.multi_augment())
+    ++e;
+  BOOST_CHECK(e > 0);
+  BOOST_CHECK(e <= 2 * sqrt(2 * (n1 + n2)));
+  Graph_matching m2 = m1;
+  BOOST_CHECK(!m2.multi_augment());
+  m2.set_r(upper_bound);
+  e = 0;
+  while (m2.multi_augment())
+    ++e;
+  BOOST_CHECK(e <= 2 * sqrt(2 * (n1 + n2)));
+  BOOST_CHECK(m2.perfect());
+  BOOST_CHECK(!m1.perfect());
 }
 
-BOOST_AUTO_TEST_CASE(global){
-    std::uniform_real_distribution<double> unif1(0.,upper_bound);
-    std::uniform_real_distribution<double> unif2(upper_bound/10000.,upper_bound/100.);
-    std::default_random_engine re;
-    std::vector< std::pair<double, double> > v1, v2;
-    for (int i = 0; i < n1; i++) {
-        double a = unif1(re);
-        double b = unif1(re);
-        double x = unif2(re);
-        double y = unif2(re);
-        v1.emplace_back(std::min(a,b), std::max(a,b));
-        v2.emplace_back(std::min(a,b)+std::min(x,y), std::max(a,b)+std::max(x,y));
-        if(i%5==0)
-            v1.emplace_back(std::min(a,b),std::min(a,b)+x);
-        if(i%3==0)
-            v2.emplace_back(std::max(a,b),std::max(a,b)+y);
-    }
-    BOOST_CHECK(bottleneck_distance(v1, v2, 0.) <= upper_bound/100.);
-    BOOST_CHECK(bottleneck_distance(v1, v2, upper_bound/10000.) <= upper_bound/100. + upper_bound/10000.);
-    BOOST_CHECK(std::abs(bottleneck_distance(v1, v2, 0.) - bottleneck_distance(v1, v2, upper_bound/10000.)) <=  upper_bound/10000.);
+BOOST_AUTO_TEST_CASE(global) {
+  std::uniform_real_distribution<double> unif1(0., upper_bound);
+  std::uniform_real_distribution<double> unif2(upper_bound / 10000., upper_bound / 100.);
+  std::default_random_engine re;
+  std::vector< std::pair<double, double> > v1, v2;
+  for (int i = 0; i < n1; i++) {
+    double a = unif1(re);
+    double b = unif1(re);
+    double x = unif2(re);
+    double y = unif2(re);
+    v1.emplace_back(std::min(a, b), std::max(a, b));
+    v2.emplace_back(std::min(a, b) + std::min(x, y), std::max(a, b) + std::max(x, y));
+    if (i % 5 == 0)
+      v1.emplace_back(std::min(a, b), std::min(a, b) + x);
+    if (i % 3 == 0)
+      v2.emplace_back(std::max(a, b), std::max(a, b) + y);
+  }
+  BOOST_CHECK(bottleneck_distance(v1, v2, 0.) <= upper_bound / 100.);
+  BOOST_CHECK(bottleneck_distance(v1, v2, upper_bound / 10000.) <= upper_bound / 100. + upper_bound / 10000.);
+  BOOST_CHECK(std::abs(bottleneck_distance(v1, v2, 0.) - bottleneck_distance(v1, v2, upper_bound / 10000.)) <= upper_bound / 10000.);
 }