Add of Bottleneck functionnality

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

Former-commit-id: 18db594b00f5065caf703e0ef026c707c985be79

14 files changed, 767 insertions, 14 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index d5dcb150..89440490 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -70,6 +70,7 @@ else()
 
   include_directories(src/common/include/)
   include_directories(src/Alpha_shapes/include/)
+  include_directories(src/Bottleneck/include/)
   include_directories(src/Contraction/include/)
   include_directories(src/Hasse_complex/include/)
   include_directories(src/Persistent_cohomology/include/)
@@ -86,6 +87,8 @@ else()
   add_subdirectory(src/Hasse_complex/example)
   add_subdirectory(src/Alpha_shapes/example)
   add_subdirectory(src/Alpha_shapes/test)
+  add_subdirectory(src/Bottleneck/example)
+  add_subdirectory(src/Bottleneck/test)
 
   # data points generator
   add_subdirectory(data/points/generator)
diff --git a/data/points/generator/CMakeLists.txt b/data/points/generator/CMakeLists.txt
index 8ff9ba52..0f2674c4 100644
--- a/data/points/generator/CMakeLists.txt
+++ b/data/points/generator/CMakeLists.txt
@@ -2,20 +2,25 @@ cmake_minimum_required(VERSION 2.6)
 project(GUDHIPointsGenerator)
 
 if(CGAL_FOUND)
-  include( ${CGAL_USE_FILE} )
+  if (NOT CGAL_VERSION VERSION_LESS 4.6.0)
+    message(STATUS "CGAL version: ${CGAL_VERSION}.")
+    include( ${CGAL_USE_FILE} )
 
-  find_package(Eigen3 3.1.0)
-  if (EIGEN3_FOUND)
-    include( ${EIGEN3_USE_FILE} )
-    include_directories (BEFORE "../../include")
+    find_package(Eigen3 3.1.0)
+    if (EIGEN3_FOUND)
+      include( ${EIGEN3_USE_FILE} )
+      include_directories (BEFORE "../../include")
 
-    add_executable ( hypergenerator hypergenerator.cpp )
-    add_test(hypergenerator_on_sphere_3000_10_5.0 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator on sphere onSphere.off 3000 10 5.0)
-    add_test(hypergenerator_on_sphere_10000_3 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator on sphere onSphere.off 10000 3)
-    add_test(hypergenerator_in_sphere_7000_12_10.8 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in sphere inSphere.off 7000 12 10.8)
-    add_test(hypergenerator_in_sphere_50000_2 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in sphere inSphere.off 50000 2)
-    # on cube is not available in CGAL
-    add_test(hypergenerator_in_cube_7000_12_10.8 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in cube inCube.off 7000 12 10.8)
-    add_test(hypergenerator_in_cube_50000_2 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in cube inCube.off 50000 3)
-  endif()
+      add_executable ( hypergenerator hypergenerator.cpp )
+      add_test(hypergenerator_on_sphere_3000_10_5.0 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator on sphere onSphere.off 3000 10 5.0)
+      add_test(hypergenerator_on_sphere_10000_3 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator on sphere onSphere.off 10000 3)
+      add_test(hypergenerator_in_sphere_7000_12_10.8 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in sphere inSphere.off 7000 12 10.8)
+      add_test(hypergenerator_in_sphere_50000_2 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in sphere inSphere.off 50000 2)
+      # on cube is not available in CGAL
+      add_test(hypergenerator_in_cube_7000_12_10.8 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in cube inCube.off 7000 12 10.8)
+      add_test(hypergenerator_in_cube_50000_2 ${CMAKE_CURRENT_BINARY_DIR}/hypergenerator in cube inCube.off 50000 3)
+      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/concept/Persistence_diagram.h b/src/Bottleneck/concept/Persistence_diagram.h
new file mode 100644
index 00000000..eaaf8bc5
--- /dev/null
+++ b/src/Bottleneck/concept/Persistence_diagram.h
@@ -0,0 +1,7 @@
+typedef typename std::pair<double,double> Diagram_point;
+
+struct Persistence_Diagram
+{
+    const_iterator<Diagram_point> cbegin() const;
+    const_iterator<Diagram_point> cend() const;
+};
diff --git a/src/Bottleneck/example/CMakeLists.txt b/src/Bottleneck/example/CMakeLists.txt
new file mode 100644
index 00000000..2ff009c4
--- /dev/null
+++ b/src/Bottleneck/example/CMakeLists.txt
@@ -0,0 +1,5 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIBottleneckExample)
+
+add_executable ( RandomDiagrams random_diagrams.cpp )
+add_test(RandomDiagrams ${CMAKE_CURRENT_BINARY_DIR}/RandomDiagrams)
diff --git a/src/Bottleneck/example/random_diagrams.cpp b/src/Bottleneck/example/random_diagrams.cpp
new file mode 100644
index 00000000..71f152a6
--- /dev/null
+++ b/src/Bottleneck/example/random_diagrams.cpp
@@ -0,0 +1,40 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "gudhi/Graph_matching.h"
+#include <iostream>
+
+using namespace Gudhi::bottleneck;
+
+int main() {
+  int n = 100;
+  std::vector< std::pair<double, double> > v1, v2;
+  for (int i = 0; i < n; i++) {
+    int a = rand() % n;
+    v1.emplace_back(a, a + rand() % (n - a));
+    int b = rand() % n;
+    v2.emplace_back(b, b + rand() % (n - b));
+  }
+  // v1 and v2 are persistence diagrams containing each 100 randoms points.
+  double b = bottleneck_distance(v1, v2, 0);
+  std::cout << b << std::endl;
+}
diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h
new file mode 100644
index 00000000..ea47e1d5
--- /dev/null
+++ b/src/Bottleneck/include/gudhi/Graph_matching.h
@@ -0,0 +1,197 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_
+
+#include <deque>
+#include <list>
+#include <vector>
+
+#include "gudhi/Layered_neighbors_finder.h"
+
+namespace Gudhi {
+
+namespace bottleneck {
+
+template<typename Persistence_diagram1, typename Persistence_diagram2>
+double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.);
+
+class Graph_matching {
+ public:
+  Graph_matching(const Persistence_diagrams_graph& g);
+  Graph_matching& operator=(const Graph_matching& m);
+  bool perfect() const;
+  bool multi_augment();
+  void set_r(double r);
+
+ private:
+  const Persistence_diagrams_graph& g;
+  double r;
+  std::vector<int> v_to_u;
+  std::list<int> unmatched_in_u;
+
+  Layered_neighbors_finder* layering() const;
+  bool augment(Layered_neighbors_finder* layered_nf, int u_start_index, int max_depth);
+  void update(std::deque<int>& path);
+};
+
+Graph_matching::Graph_matching(const Persistence_diagrams_graph& g)
+    : g(g), r(0), v_to_u(g.size()), unmatched_in_u() {
+  for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index)
+    unmatched_in_u.emplace_back(u_point_index);
+}
+
+Graph_matching& Graph_matching::operator=(const Graph_matching& m) {
+  r = m.r;
+  v_to_u = m.v_to_u;
+  unmatched_in_u = m.unmatched_in_u;
+  return *this;
+}
+
+inline bool Graph_matching::perfect() const {
+  return unmatched_in_u.empty();
+}
+
+inline bool Graph_matching::multi_augment() {
+  if (perfect())
+    return false;
+  Layered_neighbors_finder* layered_nf = layering();
+  double rn = sqrt(g.size());
+  int nblmax = layered_nf->vlayers_number()*2 + 1;
+  // verification of a necessary criterion
+  if ((unmatched_in_u.size() > rn && nblmax > rn) || nblmax == 0)
+    return false;
+  bool successful = false;
+  std::list<int>* tries = new std::list<int>(unmatched_in_u);
+  for (auto it = tries->cbegin(); it != tries->cend(); it++)
+    successful = successful || augment(layered_nf, *it, nblmax);
+  delete tries;
+  delete layered_nf;
+  return successful;
+}
+
+inline void Graph_matching::set_r(double r) {
+  this->r = r;
+}
+
+Layered_neighbors_finder* Graph_matching::layering() const {
+  bool end = false;
+  int layer = 0;
+  std::list<int> u_vertices(unmatched_in_u);
+  std::list<int> v_vertices;
+  Neighbors_finder nf(g, r);
+  Layered_neighbors_finder* layered_nf = new Layered_neighbors_finder(g, r);
+  for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index)
+    nf.add(v_point_index);
+  while (!u_vertices.empty()) {
+    for (auto it = u_vertices.cbegin(); it != u_vertices.cend(); ++it) {
+      std::list<int>* u_succ = nf.pull_all_near(*it);
+      for (auto it = u_succ->cbegin(); it != u_succ->cend(); ++it) {
+        layered_nf->add(*it, layer);
+        v_vertices.emplace_back(*it);
+      }
+      delete u_succ;
+    }
+    u_vertices.clear();
+    for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) {
+      if (v_to_u.at(*it) == null_point_index())
+        end = true;
+      else
+        u_vertices.emplace_back(v_to_u.at(*it));
+    }
+    if (end)
+      return layered_nf;
+    v_vertices.clear();
+    layer++;
+  }
+  return layered_nf;
+}
+
+bool Graph_matching::augment(Layered_neighbors_finder *layered_nf, int u_start_index, int max_depth) {
+  std::deque<int> path;
+  path.emplace_back(u_start_index);
+  // start is a point from U
+  do {
+    if (static_cast<int>(path.size()) > max_depth) {
+      path.pop_back();
+      path.pop_back();
+    }
+    if (path.empty())
+      return false;
+    int w = path.back();
+    path.emplace_back(layered_nf->pull_near(w, path.size() / 2));
+    while (path.back() == null_point_index()) {
+      path.pop_back();
+      path.pop_back();
+      if (path.empty())
+        return false;
+      path.pop_back();
+      path.emplace_back(layered_nf->pull_near(path.back(), path.size() / 2));
+    }
+    path.emplace_back(v_to_u.at(path.back()));
+  } while (path.back() != null_point_index());
+  path.pop_back();
+  update(path);
+  return true;
+}
+
+void Graph_matching::update(std::deque<int>& path) {
+  unmatched_in_u.remove(path.front());
+  for (auto it = path.cbegin(); it != path.cend(); ++it) {
+    int tmp = *it;
+    ++it;
+    v_to_u[*it] = tmp;
+  }
+}
+
+template<typename Persistence_diagram1, typename Persistence_diagram2>
+double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e) {
+  Persistence_diagrams_graph g(diag1, diag2, e);
+  std::vector<double>* sd = g.sorted_distances();
+  int idmin = 0;
+  int idmax = sd->size() - 1;
+  double alpha = pow(sd->size(), 0.25);
+  Graph_matching m(g);
+  Graph_matching biggest_unperfect = m;
+  while (idmin != idmax) {
+    int pas = static_cast<int>((idmax - idmin) / alpha);
+    m.set_r(sd->at(idmin + pas));
+    while (m.multi_augment()) {}
+    if (m.perfect()) {
+      idmax = idmin + pas;
+      m = biggest_unperfect;
+    } else {
+      biggest_unperfect = m;
+      idmin = idmin + pas + 1;
+    }
+  }
+  double b = sd->at(idmin);
+  delete sd;
+  return b;
+}
+
+}  // namespace bottleneck
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_
diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
new file mode 100644
index 00000000..de36e00b
--- /dev/null
+++ b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
@@ -0,0 +1,74 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_
+
+#include <vector>
+
+#include "Neighbors_finder.h"
+
+// Layered_neighbors_finder is a data structure used to find if a query point from U has neighbors in V in a given
+// vlayer of the vlayered persistence diagrams graph. V's points have to be added manually using their index.
+// A neighbor returned is automatically removed.
+
+namespace Gudhi {
+
+namespace bottleneck {
+
+class Layered_neighbors_finder {
+ public:
+  Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r);
+  void add(int v_point_index, int vlayer);
+  int pull_near(int u_point_index, int vlayer);
+  int vlayers_number() const;
+
+ private:
+  const Persistence_diagrams_graph& g;
+  const double r;
+  std::vector<Neighbors_finder> neighbors_finder;
+};
+
+Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r) :
+    g(g), r(r), neighbors_finder() { }
+
+inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) {
+  for (int l = neighbors_finder.size(); l <= vlayer; l++)
+    neighbors_finder.emplace_back(Neighbors_finder(g, r));
+  neighbors_finder.at(vlayer).add(v_point_index);
+}
+
+inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) {
+  if (static_cast<int> (neighbors_finder.size()) <= vlayer)
+    return null_point_index();
+  return neighbors_finder.at(vlayer).pull_near(u_point_index);
+}
+
+inline int Layered_neighbors_finder::vlayers_number() const {
+  return neighbors_finder.size();
+}
+
+}  // namespace bottleneck
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_LAYERED_NEIGHBORS_FINDER_H_
diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h
new file mode 100644
index 00000000..98256571
--- /dev/null
+++ b/src/Bottleneck/include/gudhi/Neighbors_finder.h
@@ -0,0 +1,96 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_
+
+#include <unordered_set>
+#include <list>
+
+#include "gudhi/Planar_neighbors_finder.h"
+
+namespace Gudhi {
+
+namespace bottleneck {
+
+// Neighbors_finder is a data structure used to find if a query point from U has neighbors in V
+// in the persistence diagrams graph.
+// V's points have to be added manually using their index. A neighbor returned is automatically removed.
+
+class Neighbors_finder {
+ public:
+  Neighbors_finder(const Persistence_diagrams_graph& g, double r);
+  void add(int v_point_index);
+  int pull_near(int u_point_index);
+  std::list<int>* pull_all_near(int u_point_index);
+
+ private:
+  const Persistence_diagrams_graph& g;
+  const double r;
+  Planar_neighbors_finder planar_neighbors_f;
+  std::unordered_set<int> projections_f;
+};
+
+Neighbors_finder::Neighbors_finder(const Persistence_diagrams_graph& g, double r) :
+    g(g), r(r), planar_neighbors_f(g, r), projections_f() { }
+
+inline void Neighbors_finder::add(int v_point_index) {
+  if (g.on_the_v_diagonal(v_point_index))
+    projections_f.emplace(v_point_index);
+  else
+    planar_neighbors_f.add(v_point_index);
+}
+
+inline int Neighbors_finder::pull_near(int u_point_index) {
+  int v_challenger = g.corresponding_point_in_v(u_point_index);
+  if (planar_neighbors_f.contains(v_challenger) && g.distance(u_point_index, v_challenger) < r) {
+    planar_neighbors_f.remove(v_challenger);
+    return v_challenger;
+  }
+  if (g.on_the_u_diagonal(u_point_index)) {
+    auto it = projections_f.cbegin();
+    if (it != projections_f.cend()) {
+      int tmp = *it;
+      projections_f.erase(it);
+      return tmp;
+    }
+  } else {
+    return planar_neighbors_f.pull_near(u_point_index);
+  }
+  return null_point_index();
+}
+
+inline std::list<int>* Neighbors_finder::pull_all_near(int u_point_index) {
+  std::list<int>* all_pull = planar_neighbors_f.pull_all_near(u_point_index);
+  int last_pull = pull_near(u_point_index);
+  while (last_pull != null_point_index()) {
+    all_pull->emplace_back(last_pull);
+    last_pull = pull_near(u_point_index);
+  }
+  return all_pull;
+}
+
+}  // namespace bottleneck
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_
diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
new file mode 100644
index 00000000..7e278209
--- /dev/null
+++ b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
@@ -0,0 +1,147 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_
+
+#include <vector>
+#include <set>
+#include <cmath>
+#include <utility>  // for pair<>
+#include <algorithm>  // for max
+
+namespace Gudhi {
+
+namespace bottleneck {
+
+// Diagram_point is the type of the persistence diagram's points
+typedef typename std::pair<double, double> Diagram_point;
+
+// Return the used index for encoding none of the points
+int null_point_index();
+
+// Persistence_diagrams_graph is the interface beetwen any external representation of the two persistence diagrams and
+// the bottleneck distance computation. An interface is necessary to ensure basic functions complexity.
+
+class Persistence_diagrams_graph {
+ public:
+  // Persistence_diagram1 and 2 are the types of any externals representations of persistence diagrams.
+  // They have to have an iterator over points, which have to have fields first (for birth) and second (for death).
+  template<typename Persistence_diagram1, typename Persistence_diagram2>
+  Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.);
+  Persistence_diagrams_graph();
+  bool on_the_u_diagonal(int u_point_index) const;
+  bool on_the_v_diagonal(int v_point_index) const;
+  int corresponding_point_in_u(int v_point_index) const;
+  int corresponding_point_in_v(int u_point_index) const;
+  double distance(int u_point_index, int v_point_index) const;
+  int size() const;
+  std::vector<double>* sorted_distances();
+
+ private:
+  std::vector<Diagram_point> u;
+  std::vector<Diagram_point> v;
+  Diagram_point get_u_point(int u_point_index) const;
+  Diagram_point get_v_point(int v_point_index) const;
+};
+
+inline int null_point_index() {
+  return -1;
+}
+
+template<typename Persistence_diagram1, typename Persistence_diagram2>
+Persistence_diagrams_graph::Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e)
+    : u(), v() {
+  for (auto it = diag1.cbegin(); it != diag1.cend(); ++it)
+    if (it->second - it->first > e)
+      u.emplace_back(*it);
+  for (auto it = diag2.cbegin(); it != diag2.cend(); ++it)
+    if (it->second - it->first > e)
+      v.emplace_back(*it);
+  if (u.size() < v.size())
+    swap(u, v);
+}
+
+Persistence_diagrams_graph::Persistence_diagrams_graph::Persistence_diagrams_graph()
+    : u(), v() { }
+
+inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) const {
+  return u_point_index >= static_cast<int> (u.size());
+}
+
+inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) const {
+  return v_point_index >= static_cast<int> (v.size());
+}
+
+inline int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) const {
+  return on_the_v_diagonal(v_point_index) ?
+      v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size());
+}
+
+inline int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) const {
+  return on_the_u_diagonal(u_point_index) ?
+      u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size());
+}
+
+inline double Persistence_diagrams_graph::distance(int u_point_index, int v_point_index) const {
+  // could be optimized for the case where one point is the projection of the other
+  if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index))
+    return 0;
+  Diagram_point p_u = get_u_point(u_point_index);
+  Diagram_point p_v = get_v_point(v_point_index);
+  return std::max(std::fabs(p_u.first - p_v.first), std::fabs(p_u.second - p_v.second));
+}
+
+inline int Persistence_diagrams_graph::size() const {
+  return static_cast<int> (u.size() + v.size());
+}
+
+inline std::vector<double>* Persistence_diagrams_graph::sorted_distances() {
+  // could be optimized
+  std::set<double> sorted_distances;
+  for (int u_point_index = 0; u_point_index < size(); ++u_point_index)
+    for (int v_point_index = 0; v_point_index < size(); ++v_point_index)
+      sorted_distances.emplace(distance(u_point_index, v_point_index));
+  return new std::vector<double>(sorted_distances.cbegin(), sorted_distances.cend());
+}
+
+inline Diagram_point Persistence_diagrams_graph::get_u_point(int u_point_index) const {
+  if (!on_the_u_diagonal(u_point_index))
+    return u.at(u_point_index);
+  Diagram_point projector = v.at(corresponding_point_in_v(u_point_index));
+  double x = (projector.first + projector.second) / 2;
+  return Diagram_point(x, x);
+}
+
+inline Diagram_point Persistence_diagrams_graph::get_v_point(int v_point_index) const {
+  if (!on_the_v_diagonal(v_point_index))
+    return v.at(v_point_index);
+  Diagram_point projector = u.at(corresponding_point_in_u(v_point_index));
+  double x = (projector.first + projector.second) / 2;
+  return Diagram_point(x, x);
+}
+
+}  // namespace bottleneck
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_
diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
new file mode 100644
index 00000000..4af672e4
--- /dev/null
+++ b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
@@ -0,0 +1,119 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Francois Godi
+ *
+ *    Copyright (C) 2015  INRIA Saclay (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_
+#define SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_
+
+#include <list>
+#include <iostream>
+#include <set>
+
+#include "Persistence_diagrams_graph.h"
+
+namespace Gudhi {
+
+namespace bottleneck {
+
+// Planar_neighbors_finder is a data structure used to find if a query point from U has planar neighbors in V with the
+// planar distance.
+// V's points have to be added manually using their index. A neighbor returned is automatically removed but we can also
+// remove points manually using their index.
+
+class Abstract_planar_neighbors_finder {
+ public:
+  Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r);
+  virtual ~Abstract_planar_neighbors_finder() = 0;
+  virtual void add(int v_point_index) = 0;
+  virtual void remove(int v_point_index) = 0;
+  virtual bool contains(int v_point_index) const = 0;
+  virtual int pull_near(int u_point_index) = 0;
+  virtual std::list<int>* pull_all_near(int u_point_index);
+
+ protected:
+  const Persistence_diagrams_graph& g;
+  const double r;
+};
+
+
+// Naive_pnf is a nave implementation of Abstract_planar_neighbors_finder
+
+class Naive_pnf : public Abstract_planar_neighbors_finder {
+ public:
+  Naive_pnf(const Persistence_diagrams_graph& g, double r);
+  void add(int v_point_index);
+  void remove(int v_point_index);
+  bool contains(int v_point_index) const;
+  int pull_near(int u_point_index);
+
+ private:
+  std::set<int> candidates;
+};
+
+
+// Planar_neighbors_finder is the used Abstract_planar_neighbors_finder's implementation
+typedef Naive_pnf Planar_neighbors_finder;
+
+Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r) :
+    g(g), r(r) { }
+
+inline Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() { }
+
+inline std::list<int>* Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) {
+  std::list<int>* all_pull = new std::list<int>();
+  int last_pull = pull_near(u_point_index);
+  while (last_pull != null_point_index()) {
+    all_pull->emplace_back(last_pull);
+    last_pull = pull_near(u_point_index);
+  }
+  return all_pull;
+}
+
+Naive_pnf::Naive_pnf(const Persistence_diagrams_graph& g, double r) :
+    Abstract_planar_neighbors_finder(g, r), candidates() { }
+
+inline void Naive_pnf::add(int v_point_index) {
+  candidates.emplace(v_point_index);
+}
+
+inline void Naive_pnf::remove(int v_point_index) {
+  candidates.erase(v_point_index);
+}
+
+inline bool Naive_pnf::contains(int v_point_index) const {
+  return (candidates.count(v_point_index) > 0);
+}
+
+inline int Naive_pnf::pull_near(int u_point_index) {
+  for (auto it = candidates.begin(); it != candidates.end(); ++it)
+    if (g.distance(u_point_index, *it) <= r) {
+      int tmp = *it;
+      candidates.erase(it);
+      return tmp;
+    }
+  return null_point_index();
+}
+
+}  // namespace bottleneck
+
+}  // namespace Gudhi
+
+#endif  // SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_
diff --git a/src/Bottleneck/test/CMakeLists.txt b/src/Bottleneck/test/CMakeLists.txt
new file mode 100644
index 00000000..7044372e
--- /dev/null
+++ b/src/Bottleneck/test/CMakeLists.txt
@@ -0,0 +1,21 @@
+cmake_minimum_required(VERSION 2.6)
+project(GUDHIBottleneckUnitTest)
+
+if(NOT MSVC)
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --coverage")
+  set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} --coverage")	
+  set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} --coverage")
+endif()
+
+add_executable ( BottleneckUnitTest bottleneck_unit_test.cpp )
+target_link_libraries(BottleneckUnitTest ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+
+# Unitary tests
+add_test(BottleneckUnitTest ${CMAKE_CURRENT_BINARY_DIR}/BottleneckUnitTest)
+
+if (LCOV_PATH)
+  # Lcov code coverage of unitary test
+  add_test(src/Bottleneck/lcov/coverage.log ${CMAKE_SOURCE_DIR}/scripts/check_code_coverage.sh ${CMAKE_SOURCE_DIR}/src/Bottleneck)
+endif()
+
+cpplint_add_tests("${CMAKE_SOURCE_DIR}/src/Bottleneck/include/gudhi")
diff --git a/src/Bottleneck/test/README b/src/Bottleneck/test/README
new file mode 100644
index 00000000..0e7b8673
--- /dev/null
+++ b/src/Bottleneck/test/README
@@ -0,0 +1,12 @@
+To compile:
+***********
+
+cmake .
+make
+
+To launch with details:
+***********************
+
+./BottleneckUnitTest --report_level=detailed --log_level=all 
+
+   ==> echo $? returns 0 in case of success (non-zero otherwise)
diff --git a/src/Bottleneck/test/bottleneck_unit_test.cpp b/src/Bottleneck/test/bottleneck_unit_test.cpp
new file mode 100644
index 00000000..068b8690
--- /dev/null
+++ b/src/Bottleneck/test/bottleneck_unit_test.cpp
@@ -0,0 +1,26 @@
+#define BOOST_TEST_MODULE bottleneck test
+
+#include <boost/test/included/unit_test.hpp>
+
+#include "gudhi/Graph_matching.h"
+#include <iostream>
+
+using namespace Gudhi::bottleneck;
+
+BOOST_AUTO_TEST_CASE(random_diagrams) {
+  int n = 100;
+  // Random construction
+  std::vector< std::pair<double, double> > v1, v2;
+  for (int i = 0; i < n; i++) {
+    int a = rand() % n;
+    v1.emplace_back(a, a + rand() % (n - a));
+    int b = rand() % n;
+    v2.emplace_back(b, b + rand() % (n - b));
+  }
+  // v1 and v2 are persistence diagrams containing each 100 randoms points.
+  double b = bottleneck_distance(v1, v2, 0);
+  //
+  std::cout << b << std::endl;
+  const double EXPECTED_DISTANCE = 98.5;
+  BOOST_CHECK(b == EXPECTED_DISTANCE);
+}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index d20d6d8a..ceb993fa 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -38,5 +38,6 @@ else()
   add_subdirectory(example/Contraction)                   
   add_subdirectory(example/Hasse_complex)
   add_subdirectory(example/Alpha_shapes)
+  add_subdirectory(example/Bottleneck)
 
 endif()