Merge last trunk modifications and fix new rips interface

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

Former-commit-id: 1e836cba1a372f609e94b66b6db512817a3d875e

18 files changed, 1475 insertions, 0 deletions

diff --git a/src/Bottleneck_distance/benchmark/CMakeLists.txt b/src/Bottleneck_distance/benchmark/CMakeLists.txt
new file mode 100644
index 00000000..170081ce
--- /dev/null
+++ b/src/Bottleneck_distance/benchmark/CMakeLists.txt
@@ -0,0 +1,9 @@
+cmake_minimum_required(VERSION 2.6)
+project(Bottleneck_distance_benchmark)
+
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
+  add_executable ( bottleneck_chrono bottleneck_chrono.cpp )
+  if (TBB_FOUND)
+    target_link_libraries(bottleneck_chrono ${TBB_LIBRARIES})
+  endif(TBB_FOUND)
+endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
diff --git a/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp
new file mode 100644
index 00000000..456c570b
--- /dev/null
+++ b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp
@@ -0,0 +1,62 @@
+/*    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:       Francois Godi
+ *
+ *    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
+ *    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/Bottleneck.h>
+#include <chrono>
+#include <fstream>
+#include <random>
+
+using namespace Gudhi::persistence_diagram;
+
+
+double upper_bound = 400.;  // any real > 0
+
+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.);
+    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_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));
+    result_file << n << ";" << duration.count() << ";" << b << std::endl;
+  }
+  result_file.close();
+}
diff --git a/src/Bottleneck_distance/concept/Persistence_diagram.h b/src/Bottleneck_distance/concept/Persistence_diagram.h
new file mode 100644
index 00000000..b157f22a
--- /dev/null
+++ b/src/Bottleneck_distance/concept/Persistence_diagram.h
@@ -0,0 +1,50 @@
+/*    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:       François Godi
+ *
+ *    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
+ *    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 CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_
+#define CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_
+
+namespace Gudhi {
+
+namespace persistence_diagram {
+
+/** \brief Concept of point in a persistence diagram. std::get<0>(point) must return the birth of the corresponding component and std::get<1>(point) its death.
+ * Both should be convertible to `double`.
+ * A valid implementation of this concept is std::pair<double,double>.
+ * Death should be larger than birth, death can be std::numeric_limits<double>::infinity() for components which stay alive.
+ *
+ * \ingroup bottleneck_distance
+ */
+struct DiagramPoint{};
+
+/** \brief Concept of persistence diagram. It is a range of `DiagramPoint`.
+ * std::begin(diagram) and std::end(diagram) must return corresponding iterators.
+ *
+ * \ingroup bottleneck_distance
+ */
+struct PersistenceDiagram{};
+
+}  // namespace persistence_diagram
+
+}  // namespace Gudhi
+
+#endif  // CONCEPT_BOTTLENECK_DISTANCE_PERSISTENCE_DIAGRAM_H_
diff --git a/src/Bottleneck_distance/doc/COPYRIGHT b/src/Bottleneck_distance/doc/COPYRIGHT
new file mode 100644
index 00000000..179740a6
--- /dev/null
+++ b/src/Bottleneck_distance/doc/COPYRIGHT
@@ -0,0 +1,19 @@
+The files of this directory are part of the Gudhi Library. The Gudhi library
+(Geometric Understanding in Higher Dimensions) is a generic C++ library for
+computational topology.
+
+Author(s):       François Godi
+
+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 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/>.
diff --git a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h
new file mode 100644
index 00000000..3998fe8d
--- /dev/null
+++ b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.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:       François Godi
+ *
+ *    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
+ *    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 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 {
+// needs namespace for Doxygen to link on classes
+namespace persistence_diagram {
+
+/**  \defgroup bottleneck_distance Bottleneck distance
+ * 
+ * \author    Fran&ccedil;ois Godi
+ * @{
+ * 
+ * \section bottleneckdefinition Definition
+ * 
+ * The bottleneck distance measures the similarity between two persistence diagrams. It is the shortest distance b for
+ * which there exists a perfect matching between the points of the two diagrams (completed with all the points on the
+ * diagonal in order to ignore cardinality mismatchs) such that any couple of matched points are at distance at most b.
+ *
+ * \image html perturb_pd.png On this picture, the red edges represent the matching. The bottleneck distance is the length of the longest edge.
+ *
+ */
+/** @} */  // end defgroup bottleneck_distance
+
+}  // namespace persistence_diagram
+
+}  // namespace Gudhi
+
+#endif  // DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_
diff --git a/src/Bottleneck_distance/doc/perturb_pd.png b/src/Bottleneck_distance/doc/perturb_pd.png
new file mode 100644
index 00000000..be638de0
--- /dev/null
+++ b/src/Bottleneck_distance/doc/perturb_pd.png
diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt
new file mode 100644
index 00000000..b0a19f8b
--- /dev/null
+++ b/src/Bottleneck_distance/example/CMakeLists.txt
@@ -0,0 +1,20 @@
+cmake_minimum_required(VERSION 2.6)
+project(Bottleneck_distance_examples)
+
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
+  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)
+
+  add_executable (alpha_rips_persistence_bottleneck_distance alpha_rips_persistence_bottleneck_distance.cpp)
+  target_link_libraries(alpha_rips_persistence_bottleneck_distance ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+  if (TBB_FOUND)
+    target_link_libraries(bottleneck_read_file_example ${TBB_LIBRARIES})
+    target_link_libraries(bottleneck_basic_example ${TBB_LIBRARIES})
+    target_link_libraries(alpha_rips_persistence_bottleneck_distance ${TBB_LIBRARIES})
+  endif(TBB_FOUND)
+  add_test(alpha_rips_persistence_bottleneck_distance ${CMAKE_CURRENT_BINARY_DIR}/alpha_rips_persistence_bottleneck_distance
+      ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off -r 0.15 -m 0.12 -d 3 -p 3)
+
+endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
diff --git a/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp b/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp
new file mode 100644
index 00000000..fd9f0858
--- /dev/null
+++ b/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp
@@ -0,0 +1,190 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2017  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
+ *    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/Alpha_complex.h>
+#include <gudhi/Rips_complex.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Points_off_io.h>
+#include <gudhi/Bottleneck.h>
+
+#include <CGAL/Epick_d.h>
+
+#include <boost/program_options.hpp>
+
+#include <string>
+#include <vector>
+#include <limits>  // infinity
+#include <utility>  // for pair
+#include <algorithm>  // for transform
+
+
+// Types definition
+using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = Simplex_tree::Filtration_value;
+using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
+using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
+using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >;
+using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >;
+using Point_d = Kernel::Point_d;
+using Points_off_reader = Gudhi::Points_off_reader<Point_d>;
+
+void program_options(int argc, char * argv[]
+                     , std::string & off_file_points
+                     , Filtration_value & threshold
+                     , int & dim_max
+                     , int & p
+                     , Filtration_value & min_persistence);
+
+static inline std::pair<double, double> compute_root_square(std::pair<double, double> input) {
+  return std::make_pair(std::sqrt(input.first), std::sqrt(input.second));
+}
+
+int main(int argc, char * argv[]) {
+  std::string off_file_points;
+  Filtration_value threshold;
+  int dim_max;
+  int p;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, off_file_points, threshold, dim_max, p, min_persistence);
+
+  Points_off_reader off_reader(off_file_points);
+
+  // --------------------------------------------
+  // Rips persistence
+  // --------------------------------------------
+  Rips_complex rips_complex(off_reader.get_point_cloud(), threshold, Euclidean_distance());
+
+  // Construct the Rips complex in a Simplex Tree
+  Simplex_tree rips_stree;
+
+  rips_complex.create_complex(rips_stree, dim_max);
+  std::cout << "The Rips complex contains " << rips_stree.num_simplices() << " simplices and has dimension "
+            << rips_stree.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  rips_stree.initialize_filtration();
+
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology rips_pcoh(rips_stree);
+  // initializes the coefficient field for homology
+  rips_pcoh.init_coefficients(p);
+  rips_pcoh.compute_persistent_cohomology(min_persistence);
+
+  // rips_pcoh.output_diagram();
+
+  // --------------------------------------------
+  // Alpha persistence
+  // --------------------------------------------
+  Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex(off_reader.get_point_cloud());
+
+  Simplex_tree alpha_stree;
+  alpha_complex.create_complex(alpha_stree, threshold * threshold);
+  std::cout << "The Alpha complex contains " << alpha_stree.num_simplices() << " simplices and has dimension "
+            << alpha_stree.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  alpha_stree.initialize_filtration();
+
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology alpha_pcoh(alpha_stree);
+  // initializes the coefficient field for homology
+  alpha_pcoh.init_coefficients(p);
+  alpha_pcoh.compute_persistent_cohomology(min_persistence * min_persistence);
+
+  // alpha_pcoh.output_diagram();
+
+  // --------------------------------------------
+  // Bottleneck distance between both persistence
+  // --------------------------------------------
+  double max_b_distance {};
+  for (int dim = 0; dim < dim_max; dim ++) {
+    std::vector< std::pair< Filtration_value , Filtration_value > > rips_intervals;
+    std::vector< std::pair< Filtration_value , Filtration_value > > alpha_intervals;
+    rips_intervals = rips_pcoh.intervals_in_dimension(dim);
+    alpha_intervals = alpha_pcoh.intervals_in_dimension(dim);
+    std::transform(alpha_intervals.begin(), alpha_intervals.end(), alpha_intervals.begin(), compute_root_square);
+
+    double bottleneck_distance = Gudhi::persistence_diagram::bottleneck_distance(rips_intervals, alpha_intervals);
+    std::cout << "In dimension " << dim << ", bottleneck distance = " << bottleneck_distance << std::endl;
+    if (bottleneck_distance > max_b_distance)
+      max_b_distance = bottleneck_distance;
+  }
+  std::cout << "================================================================================" << std::endl;
+  std::cout << "Bottleneck distance is " << max_b_distance << std::endl;
+
+  return 0;
+}
+
+void program_options(int argc, char * argv[]
+                     , std::string & off_file_points
+                     , Filtration_value & threshold
+                     , int & dim_max
+                     , int & p
+                     , Filtration_value & min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()
+      ("input-file", po::value<std::string>(&off_file_points),
+       "Name of an OFF file containing a point set.\n");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()
+      ("help,h", "produce help message")
+      ("max-edge-length,r",
+       po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()),
+       "Maximal length of an edge for the Rips complex construction.")
+      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
+       "Maximal dimension of the Rips complex we want to compute.")
+      ("field-charac,p", po::value<int>(&p)->default_value(11),
+       "Characteristic p of the coefficient field Z/pZ for computing homology.")
+      ("min-persistence,m", po::value<Filtration_value>(&min_persistence),
+       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).
+            options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::cout << "of a Rips complex defined on a set of input points.\n \n";
+    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
+    std::cout << "   p   dim b d \n";
+    std::cout << "where dim is the dimension of the homological feature,\n";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp
new file mode 100644
index 00000000..d0ca4e20
--- /dev/null
+++ b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp
@@ -0,0 +1,50 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Authors:       Francois Godi, small modifications by Pawel Dlotko
+ *
+ *    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
+ *    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/Bottleneck.h>
+
+#include <iostream>
+#include <vector>
+#include <utility>  // for pair
+#include <limits>  // for numeric_limits
+
+int main() {
+  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());
+
+  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);
+
+  std::cout << "Bottleneck distance = " << b << std::endl;
+
+  b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1);
+
+  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
new file mode 100644
index 00000000..bde05825
--- /dev/null
+++ b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp
@@ -0,0 +1,72 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Authors:       Francois Godi, small modifications by Pawel Dlotko
+ *
+ *    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
+ *    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 CGAL_HAS_THREADS
+
+#include <gudhi/Bottleneck.h>
+#include <iostream>
+#include <vector>
+#include <utility>  // for pair
+#include <fstream>
+#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::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
+
+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]);
+  double tolerance = 0.;
+  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
new file mode 100644
index 00000000..b90a0ee0
--- /dev/null
+++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h
@@ -0,0 +1,115 @@
+/*    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:       Francois Godi
+ *
+ *    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
+ *    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 BOTTLENECK_H_
+#define BOTTLENECK_H_
+
+#include <gudhi/Graph_matching.h>
+
+#include <vector>
+#include <algorithm>  // for max
+#include <limits>  // for numeric_limits
+
+#include <cmath>
+
+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;
+    }
+  }
+  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;
+    }
+  }
+  return sd.at(lower_bound_i);
+}
+
+/** \brief Function to compute the Bottleneck distance between two persistence diagrams.
+ *
+ * \tparam Persistence_diagram1,Persistence_diagram2
+ * models of the concept `PersistenceDiagram`.
+ * \param[in] e
+ * \parblock
+ * If `e` is 0, this uses an expensive algorithm to compute the exact distance.
+ *
+ * If `e` is not 0, it asks for an additive `e`-approximation, and currently
+ * also allows a small multiplicative error (the last 2 or 3 bits of the
+ * mantissa may be wrong). This version of the algorithm takes advantage of the
+ * limited precision of `double` and is usually a lot faster to compute,
+ * whatever the value of `e`.
+ *
+ * Thus, by default, `e` is the smallest positive double.
+ * \endparblock
+ *
+ * \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));
+}
+
+}  // namespace persistence_diagram
+
+}  // namespace Gudhi
+
+#endif  // BOTTLENECK_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..e1708c5b
--- /dev/null
+++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h
@@ -0,0 +1,182 @@
+/*    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:       Francois Godi
+ *
+ *    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
+ *    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 GRAPH_MATCHING_H_
+#define GRAPH_MATCHING_H_
+
+#include <gudhi/Neighbors_finder.h>
+
+#include <vector>
+#include <list>
+
+namespace Gudhi {
+
+namespace persistence_diagram {
+
+/** \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(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);
+}
+
+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;
+}
+
+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::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();
+  }
+  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;
+  }
+}
+
+
+}  // namespace persistence_diagram
+
+}  // namespace Gudhi
+
+#endif  // 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..0b2d26fe
--- /dev/null
+++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h
@@ -0,0 +1,91 @@
+/*    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:       Francois Godi
+ *
+ *    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
+ *    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 INTERNAL_POINT_H_
+#define INTERNAL_POINT_H_
+
+namespace Gudhi {
+
+namespace persistence_diagram {
+
+/** \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) {
+    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;
+}
+
+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;
+  }
+};
+
+}  // 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
new file mode 100644
index 00000000..cd5486f8
--- /dev/null
+++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
@@ -0,0 +1,172 @@
+/*    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:       Francois Godi
+ *
+ *    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
+ *    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 NEIGHBORS_FINDER_H_
+#define NEIGHBORS_FINDER_H_
+
+// Inclusion order is important for CGAL patch
+#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>
+
+#include <gudhi/Persistence_graph.h>
+#include <gudhi/Internal_point.h>
+
+#include <unordered_set>
+#include <vector>
+
+namespace Gudhi {
+
+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.
+ *
+ * \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;
+};
+
+/** \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(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));
+}
+
+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;
+}
+
+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;
+}
+
+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);
+}
+
+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 persistence_diagram
+
+}  // namespace Gudhi
+
+#endif  // NEIGHBORS_FINDER_H_
diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
new file mode 100644
index 00000000..44f4b827
--- /dev/null
+++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
@@ -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:       Francois Godi
+ *
+ *    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
+ *    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 PERSISTENCE_GRAPH_H_
+#define PERSISTENCE_GRAPH_H_
+
+#include <gudhi/Internal_point.h>
+
+#ifdef GUDHI_USE_TBB
+#include <tbb/parallel_sort.h>
+#endif
+
+#include <vector>
+#include <algorithm>
+#include <limits>  // for numeric_limits
+
+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 PersistenceDiagrams (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));
+  }
+}
+
+inline bool Persistence_graph::on_the_u_diagonal(int u_point_index) const {
+  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());
+}
+
+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());
+}
+
+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());
+}
+
+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()));
+}
+
+inline int Persistence_graph::size() const {
+  return static_cast<int> (u.size() + v.size());
+}
+
+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));
+  }
+#ifdef GUDHI_USE_TBB
+  tbb::parallel_sort(distances.begin(), distances.end());
+#else
+  std::sort(distances.begin(), distances.end());
+#endif
+  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);
+}
+
+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);
+}
+
+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;
+}
+
+}  // namespace persistence_diagram
+
+}  // namespace Gudhi
+
+#endif  // PERSISTENCE_GRAPH_H_
diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt
new file mode 100644
index 00000000..3d8e1f95
--- /dev/null
+++ b/src/Bottleneck_distance/test/CMakeLists.txt
@@ -0,0 +1,25 @@
+cmake_minimum_required(VERSION 2.6)
+project(Bottleneck_distance_tests)
+
+
+if (GCOVR_PATH)
+  # for gcovr to make coverage reports - Corbera Jenkins plugin
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage")
+endif()
+if (GPROF_PATH)
+  # for gprof to make coverage reports - Jenkins
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg")
+endif()
+
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
+  add_executable ( bottleneckUT bottleneck_unit_test.cpp )
+  target_link_libraries(bottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+  if (TBB_FOUND)
+    target_link_libraries(bottleneckUT ${TBB_LIBRARIES})
+  endif(TBB_FOUND)
+
+  # 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)
+endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
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_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
new file mode 100644
index 00000000..e39613b3
--- /dev/null
+++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
@@ -0,0 +1,167 @@
+/*    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:       Francois Godi
+ *
+ *    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
+ *    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_DYN_LINK
+#define BOOST_TEST_MODULE "bottleneck distance"
+#include <boost/test/unit_test.hpp>
+
+#include <random>
+#include <gudhi/Bottleneck.h>
+
+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
+
+
+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(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);
+}
+
+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);
+}
+
+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());
+}
+
+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.);
+}