renaming

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

Former-commit-id: 53bcbf20e176a21dc41ef3c4d4295b1e9aa959b0

14 files changed, 1314 insertions, 0 deletions

diff --git a/src/Bottleneck_distance/concept/Diagram_point.h b/src/Bottleneck_distance/concept/Diagram_point.h
new file mode 100644
index 00000000..ce9c6d5b
--- /dev/null
+++ b/src/Bottleneck_distance/concept/Diagram_point.h
@@ -0,0 +1,6 @@
+
+
+struct Diagram_point{
+    double first;
+    double second;
+};
diff --git a/src/Bottleneck_distance/concept/Persistence_diagram.h b/src/Bottleneck_distance/concept/Persistence_diagram.h
new file mode 100644
index 00000000..4951af32
--- /dev/null
+++ b/src/Bottleneck_distance/concept/Persistence_diagram.h
@@ -0,0 +1,7 @@
+
+
+struct Persistence_Diagram
+{
+    const_iterator<Your_Diagram_point_model> cbegin() const;
+    const_iterator<Your_Diagram_point_model> cend() const;
+};
diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt
new file mode 100644
index 00000000..4e8766cc
--- /dev/null
+++ b/src/Bottleneck_distance/example/CMakeLists.txt
@@ -0,0 +1,36 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIBottleneckExample)
+
+# need CGAL 4.6
+# cmake -DCGAL_DIR=~/workspace/CGAL-4.6-beta1 ../../..
+if(CGAL_FOUND)
+  if (NOT CGAL_VERSION VERSION_LESS 4.6.0)
+    message(STATUS "CGAL version: ${CGAL_VERSION}.")
+
+    include( ${CGAL_USE_FILE} )
+
+    if(NOT MSVC)
+      include(CheckCXXCompilerFlag)
+      CHECK_CXX_COMPILER_FLAG(-std=c++11 COMPILER_SUPPORTS_CXX11)
+      if(COMPILER_SUPPORTS_CXX11)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+      endif()
+    endif()
+    # - End of workaround
+ 
+    find_package(Eigen3 3.1.0)
+    if (EIGEN3_FOUND)
+      message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.")
+      include( ${EIGEN3_USE_FILE} )
+      include_directories (BEFORE "../../include")
+
+      add_executable (basic_example basic.cpp)
+      #add_test(dtoffrw_tore3D ${CMAKE_CURRENT_BINARY_DIR}/dtoffrw ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off 3)
+
+    else()
+      message(WARNING "Eigen3 not found. Version 3.1.0 is required.")
+    endif()
+  else()
+    message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.")
+  endif ()
+endif()
diff --git a/src/Bottleneck_distance/example/bottleneck_example.cpp b/src/Bottleneck_distance/example/bottleneck_example.cpp
new file mode 100644
index 00000000..bbd83547
--- /dev/null
+++ b/src/Bottleneck_distance/example/bottleneck_example.cpp
@@ -0,0 +1,43 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <gudhi/Graph_matching.h>
+#include <iostream>
+
+using namespace Gudhi::Bottleneck_distance;
+
+int main() {
+    std::vector< std::pair<double,double> > v1, v2;
+
+    v1.push_back(std::pair<double,double>(2.7,3.7));
+    v1.push_back(std::pair<double,double>(9.6,14));
+    v1.push_back(std::pair<double,double>(34.2,34.974));
+
+    v2.push_back(std::pair<double,double>(2.8,4.45));
+    v2.push_back(std::pair<double,double>(9.5,14.1));
+
+
+    double b = bottleneck_distance(v1, v2);
+
+    std::cout << "Bottleneck distance = " << b << std::endl;
+
+}
diff --git a/src/Bottleneck_distance/include/CGAL/Miscellaneous.h b/src/Bottleneck_distance/include/CGAL/Miscellaneous.h
new file mode 100644
index 00000000..4de787fb
--- /dev/null
+++ b/src/Bottleneck_distance/include/CGAL/Miscellaneous.h
@@ -0,0 +1,51 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_CGAL_MISCELLANEOUS_H_
+#define SRC_BOTTLENECK_INCLUDE_CGAL_MISCELLANEOUS_H_
+
+#include <gudhi/Internal_point.h>
+
+namespace CGAL {
+
+typedef Gudhi::bipartite_graph_matching::Internal_point Internal_point;
+
+template <>
+struct Kernel_traits<Internal_point> {
+    struct Kernel {
+        typedef double FT;
+        typedef double RT;
+    };
+};
+
+
+struct Construct_coord_iterator {
+    typedef  const double* result_type;
+    const double* operator()(const Internal_point& p) const
+    { return static_cast<const double*>(p.vec); }
+    const double* operator()(const Internal_point& p, int)  const
+    { return static_cast<const double*>(p.vec+2); }
+};
+
+} //namespace CGAL
+
+#endif  // SRC_BOTTLENECK_INCLUDE_CGAL_MISCELLANEOUS_H_
diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h
new file mode 100644
index 00000000..d8860841
--- /dev/null
+++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h
@@ -0,0 +1,233 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_
+
+#include <deque>
+
+#include <gudhi/Neighbors_finder.h>
+
+namespace Gudhi {
+
+namespace Bottleneck_distance {
+
+/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams.
+ *
+ *
+ *
+ * \ingroup bottleneck_distance
+ */
+template<typename Persistence_diagram1, typename Persistence_diagram2>
+double compute(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e = 0.);
+
+/** \internal \brief Structure representing a graph matching. The graph is a Persistence_diagrams_graph.
+ *
+ * \ingroup bottleneck_distance
+ */
+class Graph_matching {
+public:
+    /** \internal \brief Constructor constructing an empty matching. */
+    explicit Graph_matching();
+    /** \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:
+    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. */
+    std::shared_ptr<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::deque<int> & path);
+};
+
+inline Graph_matching::Graph_matching()
+    : 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);
+}
+
+inline Graph_matching& Graph_matching::operator=(const Graph_matching& m) {
+    r = m.r;
+    v_to_u = m.v_to_u;
+    unmatched_in_u = m.unmatched_in_u;
+    return *this;
+}
+
+inline bool Graph_matching::perfect() const {
+    return unmatched_in_u.empty();
+}
+
+inline bool Graph_matching::multi_augment() {
+    if (perfect())
+        return false;
+    Layered_neighbors_finder layered_nf = *layering();
+    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;
+}
+
+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::deque<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 std::shared_ptr<Layered_neighbors_finder> Graph_matching::layering() const {
+    std::list<int> u_vertices(unmatched_in_u);
+    std::list<int> v_vertices;
+    Neighbors_finder nf(r);
+    for (int v_point_index = 0; v_point_index < G::size(); ++v_point_index)
+        nf.add(v_point_index);
+    std::shared_ptr<Layered_neighbors_finder> layered_nf(new Layered_neighbors_finder(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::shared_ptr<std::list<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();
+    }
+    return layered_nf;
+}
+
+inline void Graph_matching::update(std::deque<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;
+    }
+}
+
+template<typename Persistence_diagram1, typename Persistence_diagram2>
+double compute_exactly(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2) {
+    G::initialize(diag1, diag2, 0.);
+    std::shared_ptr< std::vector<double> > sd(G::sorted_distances());
+    int idmin = 0;
+    int idmax = sd->size() - 1;
+    // alpha can be modified, this will change the complexity
+    double alpha = pow(sd->size(), 0.25);
+    Graph_matching m;
+    Graph_matching biggest_unperfect;
+    while (idmin != idmax) {
+        int step = static_cast<int>((idmax - idmin) / alpha);
+        m.set_r(sd->at(idmin + step));
+        while (m.multi_augment());
+        //The above while compute a maximum matching (according to the r setted before)
+        if (m.perfect()) {
+            idmax = idmin + step;
+            m = biggest_unperfect;
+        } else {
+            biggest_unperfect = m;
+            idmin = idmin + step + 1;
+        }
+    }
+    return sd->at(idmin);
+}
+
+template<typename Persistence_diagram1, typename Persistence_diagram2>
+double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) {
+    if(e< std::numeric_limits<double>::min())
+        return compute_exactly(diag1, diag2);
+    G::initialize(diag1, diag2, e);
+    double d = 0.;
+    double f = G::diameter();
+    while (f-d > e){
+        Graph_matching m;
+        m.set_r((d+f)/2.);
+        while (m.multi_augment());
+        //The above while compute a maximum matching (according to the r setted before)
+        if (m.perfect()) 
+            f = (d+f)/2.;
+         else 
+            d= (d+f)/2.;
+    }
+    return d;
+}
+
+}  // namespace Bottleneck_distance
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_
diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h
new file mode 100644
index 00000000..2080900d
--- /dev/null
+++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h
@@ -0,0 +1,74 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_INTERNAL_POINT_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_INTERNAL_POINT_H_
+
+namespace Gudhi {
+
+namespace Bottleneck_distance {
+
+/** \internal \brief Returns the used index for encoding none of the points */
+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 = null_point_index()) { 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); }
+/* 
+Useless
+    friend void swap(Internal_point& a, Internal_point& b)
+    {
+        using std::swap;
+        double x_tmp = a.vec[0];
+        double y_tmp = a.vec[1];
+        int pi_tmp = a.point_index;
+        a.vec[0] = b.vec[0];
+        a.vec[1] = b.vec[1];
+        a.point_index = b.point_index;
+        b.vec[0] =  x_tmp;
+        b.vec[1] = y_tmp;
+        b.point_index = pi_tmp;
+    }
+*/
+};
+
+inline int null_point_index() {
+    return -1;
+}
+
+}  // namespace Bottleneck_distance
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_INTERNAL_POINT_H_
diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
new file mode 100644
index 00000000..7ad79126
--- /dev/null
+++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
@@ -0,0 +1,154 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_
+
+#include <unordered_set>
+#include <gudhi/Planar_neighbors_finder.h>
+
+namespace Gudhi {
+
+namespace Bottleneck_distance {
+
+/** \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.
+ *
+ * \ingroup bottleneck_distance
+ */
+class Neighbors_finder {
+public:
+    /** \internal \brief Constructor taking the near distance definition as parameter. */
+    Neighbors_finder(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::shared_ptr< std::list<int> > pull_all_near(int u_point_index);
+
+private:
+    const double r;
+    Planar_neighbors_finder planar_neighbors_f;
+    std::unordered_set<int> projections_f;
+    void remove(int v_point_index);
+    bool contains(int v_point_index);
+};
+
+/** \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.
+ *
+ * \ingroup bottleneck_distance
+ */
+class Layered_neighbors_finder {
+public:
+    /** \internal \brief Constructor taking the near distance definition as parameter. */
+    Layered_neighbors_finder(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 double r;
+    std::vector<std::shared_ptr<Neighbors_finder>> neighbors_finder;
+};
+
+inline Neighbors_finder::Neighbors_finder(double r) :
+    r(r), planar_neighbors_f(r), projections_f() { }
+
+inline void Neighbors_finder::add(int v_point_index) {
+    if (G::on_the_v_diagonal(v_point_index))
+        projections_f.emplace(v_point_index);
+    else
+        planar_neighbors_f.add(v_point_index);
+}
+
+inline void Neighbors_finder::remove(int v_point_index) {
+    if(v_point_index == null_point_index())
+        return;
+    if (G::on_the_v_diagonal(v_point_index))
+        projections_f.erase(v_point_index);
+    else
+        planar_neighbors_f.remove(v_point_index);
+}
+
+inline bool Neighbors_finder::contains(int v_point_index) {
+    return planar_neighbors_f.contains(v_point_index) || (projections_f.count(v_point_index)>0);
+}
+
+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();
+    else if (contains(c) && (G::distance(u_point_index, c) <= r))
+        //Is the query point near to its projection ?
+        tmp = c;
+    else
+        //Is the query point near to a V point in the plane ?
+        tmp = planar_neighbors_f.pull_near(u_point_index);
+    remove(tmp);
+    return tmp;
+}
+
+inline std::shared_ptr< std::list<int> > Neighbors_finder::pull_all_near(int u_point_index) {
+    std::shared_ptr< std::list<int> > all_pull(planar_neighbors_f.pull_all_near(u_point_index));
+    int last_pull = pull_near(u_point_index);
+    while (last_pull != null_point_index()) {
+        all_pull->emplace_back(last_pull);
+        last_pull = pull_near(u_point_index);
+    }
+    return all_pull;
+}
+
+inline Layered_neighbors_finder::Layered_neighbors_finder(double r) :
+    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::shared_ptr<Neighbors_finder>(new Neighbors_finder(r)));
+    neighbors_finder.at(vlayer)->add(v_point_index);
+}
+
+inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) {
+    if (static_cast<int> (neighbors_finder.size()) <= vlayer)
+        return null_point_index();
+    return neighbors_finder.at(vlayer)->pull_near(u_point_index);
+}
+
+inline int Layered_neighbors_finder::vlayers_number() const {
+    return static_cast<int>(neighbors_finder.size());
+}
+
+}  // namespace Bottleneck_distance
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_
diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h
new file mode 100644
index 00000000..e8a46f1c
--- /dev/null
+++ b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h
@@ -0,0 +1,168 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_
+
+#include <vector>
+#include <set>
+#include <cmath>
+#include <utility>
+#include <algorithm>
+#include <math.h>
+#include <memory>
+#include <gudhi/Internal_point.h>
+
+namespace Gudhi {
+
+namespace Bottleneck_distance {
+
+
+/** \internal \brief Structure representing an euclidean bipartite graph containing
+ *  the points from the two persistence diagrams (including the projections).
+ *
+ * \ingroup bottleneck_distance
+ */
+class Persistence_diagrams_graph {
+public:
+    /** \internal \brief Initializer taking 2 Point (concept) ranges as parameters. */
+    template<typename Persistence_diagram1, typename Persistence_diagram2>
+    static void initialize(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 ? */
+    static bool on_the_u_diagonal(int u_point_index);
+    /** \internal \brief Is the given point from V the projection of a point in U ? */
+    static bool on_the_v_diagonal(int v_point_index);
+    /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */
+    static int corresponding_point_in_u(int v_point_index);
+    /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */
+    static int corresponding_point_in_v(int u_point_index);
+    /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */
+    static double distance(int u_point_index, int v_point_index);
+    /** \internal \brief Returns size = |U| = |V|. */
+    static int size();
+    /** \internal \brief Returns the O(n^2) sorted distances between the points. */
+    static std::shared_ptr< std::vector<double> > sorted_distances();
+    /** \internal \brief Returns an upper bound of the diameter of the convex hull */
+    static double diameter();
+
+private:
+    static std::vector<Internal_point> u;
+    static std::vector<Internal_point> v;
+    static Internal_point get_u_point(int u_point_index);
+    static Internal_point get_v_point(int v_point_index);
+
+    friend class Naive_pnf;
+    friend class Cgal_pnf;
+};
+
+/** \internal \typedef \brief Shorter alias */
+typedef Persistence_diagrams_graph G;
+
+// static initialization
+std::vector<Internal_point> G::u = [] {return std::vector<Internal_point>();}();
+std::vector<Internal_point> G::v = [] {return std::vector<Internal_point>();}();
+
+template<typename Persistence_diagram1, typename Persistence_diagram2>
+inline void G::initialize(const Persistence_diagram1 &diag1,
+                          const Persistence_diagram2 &diag2, double e){
+    u.clear();
+    v.clear();
+    for (auto it = diag1.cbegin(); it != diag1.cend(); ++it)
+        if (it->second - it->first > e)
+            u.push_back(Internal_point(it->first, it->second, u.size()));
+    for (auto it = diag2.cbegin(); it != diag2.cend(); ++it)
+        if (it->second - it->first > e)
+            v.push_back(Internal_point(it->first, it->second, v.size()));
+    if (u.size() < v.size())
+        swap(u, v);
+}
+
+inline bool G::on_the_u_diagonal(int u_point_index) {
+    return u_point_index >= static_cast<int> (u.size());
+}
+
+inline bool G::on_the_v_diagonal(int v_point_index) {
+    return v_point_index >= static_cast<int> (v.size());
+}
+
+inline int G::corresponding_point_in_u(int v_point_index) {
+    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 G::corresponding_point_in_v(int u_point_index) {
+    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 G::distance(int u_point_index, int v_point_index) {
+    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 G::size() {
+    return static_cast<int> (u.size() + v.size());
+}
+
+inline std::shared_ptr< std::vector<double> > G::sorted_distances() {
+    // could be optimized
+    std::set<double> sorted_distances;
+    for (int u_point_index = 0; u_point_index < size(); ++u_point_index)
+        for (int v_point_index = 0; v_point_index < size(); ++v_point_index)
+            sorted_distances.emplace(distance(u_point_index, v_point_index));
+    std::shared_ptr< std::vector<double> > sd_up(new std::vector<double>(sorted_distances.cbegin(), sorted_distances.cend()));
+    return sd_up;
+}
+
+inline Internal_point G::get_u_point(int 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 G::get_v_point(int 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 G::diameter() {
+    double max = 0.;
+    for(it = u.cbegin(); it != u.cend(); it++)
+        max = std::max(max,it->y());
+    for(it = v.cbegin(); it != v.cend(); it++)
+        max = std::max(max,it->y());
+    return max;
+}
+
+}  // namespace Bottleneck_distance
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_
diff --git a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h
new file mode 100644
index 00000000..c0b6383f
--- /dev/null
+++ b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h
@@ -0,0 +1,218 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_
+
+#include <list>
+#include <map>
+#include <CGAL/Search_traits.h>
+#include <CGAL/Orthogonal_incremental_neighbor_search.h>
+#include <CGAL/Weighted_Minkowski_distance.h>
+#include <CGAL/Miscellaneous.h>
+#include <gudhi/Persistence_diagrams_graph.h>
+
+
+namespace Gudhi {
+
+namespace Bottleneck_distance {
+
+/** \internal \brief Structure used to find any point in V near (according to the planar distance) to a query point from U.
+ *
+ * V points have to be added manually using their index and before the first remove/pull. A neighbor pulled is automatically removed. but we can also
+ * remove points manually using their index.
+ *
+ * \ingroup bottleneck_distance
+ */
+class Naive_pnf {
+public:
+    /** \internal \brief Constructor taking the near distance definition as parameter. */
+    Naive_pnf(double r_);
+    /** \internal \brief A point added will be possibly pulled. */
+    void add(int v_point_index);
+    /** \internal \brief A point manually removed will no longer be possibly pulled. */
+    void remove(int v_point_index);
+    /** \internal \brief Can the point given as parameter be returned ? */
+    bool contains(int v_point_index) const;
+    /** \internal \brief Provide 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 Provide and remove all the V points near to the U point given as parameter. */
+    virtual std::shared_ptr< std::list<int> > pull_all_near(int u_point_index);
+
+private:
+    double r;
+    std::pair<int,int> get_v_key(int v_point_index) const;
+    std::multimap<std::pair<int,int>,int> grid;
+};
+
+class Cgal_pnf {
+
+    typedef CGAL::Dimension_tag<2> D;
+    typedef CGAL::Search_traits<double, Internal_point, const double*, CGAL::Construct_coord_iterator, D> Traits;
+    typedef CGAL::Weighted_Minkowski_distance<Traits> Distance;
+    typedef CGAL::Orthogonal_incremental_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. */
+    Cgal_pnf(double r_);
+    /** \internal \brief A point added will be possibly pulled. */
+    void add(int v_point_index);
+    /** \internal \brief A point manually removed will no longer be possibly pulled. */
+    void remove(int v_point_index);
+    /** \internal \brief Can the point given as parameter be returned ? */
+    bool contains(int v_point_index) const;
+    /** \internal \brief Provide 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 Provide and remove all the V points near to the U point given as parameter. */
+    virtual std::shared_ptr< std::list<int> > pull_all_near(int u_point_index);
+
+private:
+    double r;
+    std::set<int> contents;
+    Kd_tree kd_t;
+};
+
+/** \internal \typedef \brief Planar_neighbors_finder is the used implementation. */
+typedef Naive_pnf Planar_neighbors_finder;
+
+inline Naive_pnf::Naive_pnf(double r_) :
+    r(r_), grid() { }
+
+
+inline std::pair<int,int> Naive_pnf::get_v_key(int v_point_index) const{
+    Internal_point v_point = G::get_v_point(v_point_index);
+    return std::make_pair(static_cast<int>(v_point.x()/r), static_cast<int>(v_point.y()/r));
+}
+
+inline void Naive_pnf::add(int v_point_index) {
+    grid.emplace(get_v_key(v_point_index),v_point_index);
+}
+
+inline void Naive_pnf::remove(int v_point_index) {
+    if(v_point_index != null_point_index())
+        for(auto it = grid.find(get_v_key(v_point_index)); it!=grid.end(); it++)
+            if(it->second==v_point_index){
+                grid.erase(it);
+                return;
+            }
+}
+
+inline bool Naive_pnf::contains(int v_point_index) const {
+    if(v_point_index == null_point_index())
+        return false;
+    for(auto it = grid.find(get_v_key(v_point_index)); it!=grid.end(); it++)
+        if(it->second==v_point_index)
+            return true;
+    return false;
+}
+
+inline int Naive_pnf::pull_near(int u_point_index) {
+    Internal_point u_point = G::get_u_point(u_point_index);
+    int i0 = static_cast<int>(u_point.x()/r);
+    int j0 = static_cast<int>(u_point.y()/r);
+    for(int i = 1; i<= 3; i++)
+        for(int j = 1; j<= 3; j++)
+            for(auto it = grid.find(std::make_pair(i0 +(i%3)-1, j0+(j%3)-1)); it!=grid.end(); it++)
+                if (G::distance(u_point_index, it->second) <= r) {
+                    int tmp = it->second;
+                    grid.erase(it);
+                    return tmp;
+                }
+    return null_point_index();
+}
+
+inline std::shared_ptr< std::list<int> > Naive_pnf::pull_all_near(int u_point_index) {
+    std::shared_ptr< std::list<int> > all_pull(new std::list<int>);
+    Internal_point u_point = G::get_u_point(u_point_index);
+    int i0 = static_cast<int>(u_point.x()/r);
+    int j0 = static_cast<int>(u_point.y()/r);
+    for(int i = 1; i<= 3; i++)
+        for(int j = 1; j<= 3; j++)
+            for(auto it = grid.find(std::make_pair(i0 +(i%3)-1, j0+(j%3)-1)); it!=grid.end();)
+                if (G::distance(u_point_index, it->second) <= r) {
+                    all_pull->emplace_back(it->second);
+                    it = grid.erase(it);
+                }
+                else it++;
+    return all_pull;
+}
+
+
+/** \internal \brief Constructor taking the near distance definition as parameter. */
+inline Cgal_pnf::Cgal_pnf(double r_)
+    : r(r_), contents(), kd_t() {}
+
+
+/** \internal \brief A point added will be possibly pulled. */
+inline void Cgal_pnf::add(int v_point_index){
+    if(v_point_index == null_point_index())
+        return;
+    contents.insert(v_point_index);
+    kd_t.insert(G::get_v_point(v_point_index));
+}
+
+/** \internal \brief A point manually removed will no longer be possibly pulled. */
+inline void Cgal_pnf::remove(int v_point_index){
+    if(contains(v_point_index)){
+        contents.erase(v_point_index);
+        kd_t.remove(G::get_v_point(v_point_index));
+    }
+}
+
+/** \internal \brief Can the point given as parameter be returned ? */
+inline bool Cgal_pnf::contains(int v_point_index) const{
+    if(v_point_index == null_point_index())
+        return false;
+    return contents.count(v_point_index)>0;
+}
+
+/** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */
+inline int Cgal_pnf::pull_near(int u_point_index){
+    Internal_point u_point = G::get_u_point(u_point_index);
+    std::vector<double> w = {1., 1.};
+    K_neighbor_search search(kd_t, u_point, 0., true, Distance(0, 2, w));
+    auto it = search.begin();
+    if(it==search.end() || G::distance(u_point_index, it->first.point_index) > r)
+        return null_point_index();
+    int tmp = it->first.point_index;
+    remove(tmp);
+    return tmp;
+}
+
+inline std::shared_ptr< std::list<int> > Cgal_pnf::pull_all_near(int u_point_index) {
+    std::shared_ptr< std::list<int> > all_pull(new std::list<int>);
+    int last_pull = pull_near(u_point_index);
+    while (last_pull != null_point_index()) {
+        all_pull->emplace_back(last_pull);
+        last_pull = pull_near(u_point_index);
+    }
+    return all_pull;
+}
+
+
+}  // namespace Bottleneck_distance
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_
diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt
new file mode 100644
index 00000000..ab0461cc
--- /dev/null
+++ b/src/Bottleneck_distance/test/CMakeLists.txt
@@ -0,0 +1,63 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIBottleneckUT)
+
+
+
+
+if(CGAL_FOUND)
+  if (NOT CGAL_VERSION VERSION_LESS 4.6.0)
+    message(STATUS "CGAL version: ${CGAL_VERSION}.")
+
+    include( ${CGAL_USE_FILE} )
+
+    if(NOT MSVC)
+      include(CheckCXXCompilerFlag)
+      CHECK_CXX_COMPILER_FLAG(-std=c++11 COMPILER_SUPPORTS_CXX11)
+      if(COMPILER_SUPPORTS_CXX11)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+      endif()
+    endif()
+    # - End of workaround
+ 
+    find_package(Eigen3 3.1.0)
+    if (EIGEN3_FOUND)
+      message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.")
+      include( ${EIGEN3_USE_FILE} )
+      include_directories (BEFORE "../../include")
+
+      if (GCOVR_PATH)
+  # for gcovr to make coverage reports - Corbera Jenkins plugin
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage")
+  set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fprofile-arcs -ftest-coverage")
+  set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fprofile-arcs -ftest-coverage")
+endif()
+if (GPROF_PATH)
+  # for gprof to make coverage reports - Jenkins
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg")
+  set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pg")
+  set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -pg")
+endif()
+  message("CMAKE_CXX_FLAGS = ${CMAKE_CXX_FLAGS}")
+  message("CMAKE_CXX_FLAGS_DEBUG = ${CMAKE_CXX_FLAGS_DEBUG}")	
+  message("CMAKE_CXX_FLAGS_RELEASE = ${CMAKE_CXX_FLAGS_RELEASE}")
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -std=c++11 -Wall -Wpedantic -Wsign-compare")
+  set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -ggdb -O0")	
+  set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}")
+
+add_executable ( BottleneckUT bottleneck_unit_test.cpp )
+target_link_libraries(BottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+
+# Unitary tests
+add_test(NAME BottleneckUT 
+  COMMAND ${CMAKE_CURRENT_BINARY_DIR}/BottleneckUT
+  # XML format for Jenkins xUnit plugin 
+  --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/BottleneckUT.xml --log_level=test_suite --report_level=no)
+
+
+    else()
+      message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.")
+    endif()
+  else()
+    message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile bipartite graphs matching feature. Version 4.6.0 is required.")
+  endif ()
+endif()
diff --git a/src/Bottleneck_distance/test/README b/src/Bottleneck_distance/test/README
new file mode 100644
index 00000000..0e7b8673
--- /dev/null
+++ b/src/Bottleneck_distance/test/README
@@ -0,0 +1,12 @@
+To compile:
+***********
+
+cmake .
+make
+
+To launch with details:
+***********************
+
+./BottleneckUnitTest --report_level=detailed --log_level=all 
+
+   ==> echo $? returns 0 in case of success (non-zero otherwise)
diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp
new file mode 100644
index 00000000..8f187a91
--- /dev/null
+++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp
@@ -0,0 +1,61 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <gudhi/Graph_matching.h>
+#include <chrono>
+#include <fstream>
+
+using namespace Gudhi::Bottleneck_distance;
+
+
+double upper_bound = 400.; // any real >0
+
+int main(){
+    std::ofstream objetfichier;
+    objetfichier.open("results.csv", std::ios::out);
+
+    for(int n=0; n<=2500; n+=250){
+        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);
+        }
+        std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
+        double b = bottleneck_approx(v1,v2, 0.);
+        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();
+}
diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
new file mode 100644
index 00000000..33ac3f1b
--- /dev/null
+++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
@@ -0,0 +1,188 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+#define BOOST_TEST_MODULE bottleneck test
+
+#include <boost/test/included/unit_test.hpp>
+#include <random>
+#include <gudhi/Graph_matching.h>
+
+using namespace Gudhi::Bottleneck_distance;
+
+int n1 = 81; // a natural number >0
+int n2 = 180; // a natural number >0
+double upper_bound = 406.43; // any real >0
+
+BOOST_AUTO_TEST_CASE(persistence_diagrams_graph){
+    // Random construction
+    std::uniform_real_distribution<double> unif(0.,upper_bound);
+    std::default_random_engine re;
+    std::vector< std::pair<double, double> > v1, v2;
+    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));
+    }
+    G::initialize(v1, v2, 0.);
+    std::shared_ptr< 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))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n1-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n2-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n2))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,(n1+n2)-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,0))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n1-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n2-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n2))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,(n1+n2)-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,0))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n1-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n2-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n2))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,(n1+n2)-1))==1);
+}
+
+BOOST_AUTO_TEST_CASE(planar_neighbors_finder) {
+    Planar_neighbors_finder pnf(1.);
+    for(int v_point_index=0; v_point_index<n1; v_point_index+=2)
+        pnf.add(v_point_index);
+    //
+    BOOST_CHECK(pnf.contains(0));
+    BOOST_CHECK(!pnf.contains(1));
+    BOOST_CHECK(pnf.contains(2));
+    BOOST_CHECK(!pnf.contains(3));
+    //
+    pnf.remove(0);
+    pnf.remove(1);
+    //
+    BOOST_CHECK(!pnf.contains(0));
+    BOOST_CHECK(!pnf.contains(1));
+    BOOST_CHECK(pnf.contains(2));
+    BOOST_CHECK(!pnf.contains(3));
+    //
+    int v_point_index_1 = pnf.pull_near(n2/2);
+    BOOST_CHECK((v_point_index_1 == -1) || ((G::distance(n2/2,v_point_index_1)<=1.)));
+    BOOST_CHECK(!pnf.contains(v_point_index_1));
+    std::list<int> l = *pnf.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 = pnf.pull_near(n2/2);
+    BOOST_CHECK(v_point_index_2 == -1);
+}
+
+BOOST_AUTO_TEST_CASE(neighbors_finder) {
+    Neighbors_finder nf(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::list<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) {
+    Layered_neighbors_finder lnf(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) {
+    Graph_matching m1;
+    m1.set_r(0.);
+    int e  = 0;
+    while (m1.multi_augment())
+        ++e;
+    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/1000.,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) <= upper_bound/100.);
+}