Reimplement the bottleneck python wrapper with pybind11

6 files changed, 125 insertions, 97 deletions

diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index d712e189..976a8b52 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -34,6 +34,7 @@ endfunction( add_gudhi_debug_info )
 if(PYTHONINTERP_FOUND)
   if(PYBIND11_FOUND)
     add_gudhi_debug_info("Pybind11 version ${PYBIND11_VERSION}")
+    set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'bottleneck', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'hera', ")
   endif()
   if(CYTHON_FOUND)
@@ -46,7 +47,6 @@ if(PYTHONINTERP_FOUND)
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'reader_utils', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'witness_complex', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'strong_witness_complex', ")
-    set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'bottleneck', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'nerve_gic', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'subsampling', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'tangential_complex', ")
@@ -120,24 +120,25 @@ if(PYTHONINTERP_FOUND)
       set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_EIGEN3_ENABLED', ")
     endif (EIGEN3_FOUND)
 
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'off_reader', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'simplex_tree', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'rips_complex', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'cubical_complex', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'periodic_cubical_complex', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'reader_utils', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'witness_complex', ")
-    set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'strong_witness_complex', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_reader', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'simplex_tree', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'rips_complex', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'cubical_complex', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'periodic_cubical_complex', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'reader_utils', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'witness_complex', ")
+    set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'strong_witness_complex', ")
+    set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera', ")
     if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'bottleneck', ")
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'nerve_gic', ")
+      set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'bottleneck', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'nerve_gic', ")
     endif ()
     if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'alpha_complex', ")
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'subsampling', ")
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'tangential_complex', ")
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'euclidean_witness_complex', ")
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'euclidean_strong_witness_complex', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'subsampling', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'tangential_complex', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_witness_complex', ")
+      set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_strong_witness_complex', ")
     endif ()
 
     if(CGAL_FOUND)
diff --git a/src/python/gudhi/bottleneck.cc b/src/python/gudhi/bottleneck.cc
new file mode 100644
index 00000000..577e5e0b
--- /dev/null
+++ b/src/python/gudhi/bottleneck.cc
@@ -0,0 +1,51 @@
+/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ *    Author(s):       Marc Glisse
+ *
+ *    Copyright (C) 2020 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#include <gudhi/Bottleneck.h>
+
+#include <pybind11_diagram_utils.h>
+
+double bottleneck(Dgm d1, Dgm d2, double epsilon)
+{
+  // I *think* the call to request() has to be before releasing the GIL.
+  auto diag1 = numpy_to_range_of_pairs(d1);
+  auto diag2 = numpy_to_range_of_pairs(d2);
+
+  py::gil_scoped_release release;
+
+  return Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, epsilon);
+}
+
+PYBIND11_MODULE(bottleneck, m) {
+      m.attr("__license__") = "GPL v3";
+      m.def("bottleneck_distance", &bottleneck,
+          py::arg("diagram_1"), py::arg("diagram_2"),
+          py::arg("e") = (std::numeric_limits<double>::min)(),
+          R"pbdoc(
+    This function returns the point corresponding to a given vertex.
+
+    :param diagram_1: The first diagram.
+    :type diagram_1: vector[pair[double, double]]
+    :param diagram_2: The second diagram.
+    :type diagram_2: vector[pair[double, double]]
+    :param e: 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.
+    :type e: float
+    :rtype: float
+    :returns: the bottleneck distance.
+    )pbdoc");
+}
diff --git a/src/python/gudhi/bottleneck.pyx b/src/python/gudhi/bottleneck.pyx
deleted file mode 100644
index af011e88..00000000
--- a/src/python/gudhi/bottleneck.pyx
+++ /dev/null
@@ -1,48 +0,0 @@
-# This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
-# See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
-# Author(s):       Vincent Rouvreau
-#
-# Copyright (C) 2016 Inria
-#
-# Modification(s):
-#   - YYYY/MM Author: Description of the modification
-
-from cython cimport numeric
-from libcpp.vector cimport vector
-from libcpp.utility cimport pair
-import os
-
-__author__ = "Vincent Rouvreau"
-__copyright__ = "Copyright (C) 2016 Inria"
-__license__ = "GPL v3"
-
-cdef extern from "Bottleneck_distance_interface.h" namespace "Gudhi::persistence_diagram":
-    double bottleneck(vector[pair[double, double]], vector[pair[double, double]], double)
-    double bottleneck(vector[pair[double, double]], vector[pair[double, double]])
-
-def bottleneck_distance(diagram_1, diagram_2, e=None):
-    """This function returns the point corresponding to a given vertex.
-
-    :param diagram_1: The first diagram.
-    :type diagram_1: vector[pair[double, double]]
-    :param diagram_2: The second diagram.
-    :type diagram_2: vector[pair[double, double]]
-    :param e: 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.
-    :type e: float
-    :rtype: float
-    :returns: the bottleneck distance.
-    """
-    if e is None:
-        # Default value is the smallest double value (not 0, 0 is for exact version)
-        return bottleneck(diagram_1, diagram_2)
-    else:
-        # Can be 0 for exact version
-        return bottleneck(diagram_1, diagram_2, e)
diff --git a/src/python/gudhi/hera.cc b/src/python/gudhi/hera.cc
index 5aec1806..ea80a9a8 100644
--- a/src/python/gudhi/hera.cc
+++ b/src/python/gudhi/hera.cc
@@ -8,39 +8,9 @@
  *      - YYYY/MM Author: Description of the modification
  */
 
-#include <pybind11/pybind11.h>
-#include <pybind11/numpy.h>
-
-#include <boost/range/counting_range.hpp>
-#include <boost/range/adaptor/transformed.hpp>
-
 #include <wasserstein.h> // Hera
 
-#include <utility>
-
-namespace py = pybind11;
-typedef py::array_t<double> Dgm;
-
-// Get m[i,0] and m[i,1] as a pair
-static auto pairify(void* p, ssize_t h, ssize_t w) {
-  return [=](ssize_t i){
-    char* birth = (char*)p + i * h;
-    char* death = birth + w;
-    return std::make_pair(*(double*)birth, *(double*)death);
-  };
-}
-
-inline auto numpy_to_range_of_pairs(py::array_t<double> dgm) {
-  py::buffer_info buf = dgm.request();
-  // shape (n,2) or (0) for empty
-  if((buf.ndim!=2 || buf.shape[1]!=2) && (buf.ndim!=1 || buf.shape[0]!=0))
-    throw std::runtime_error("Diagram must be an array of size n x 2");
-  // In the case of shape (0), avoid reading non-existing strides[1] even if we won't use it.
-  ssize_t stride1 = buf.ndim == 2 ? buf.strides[1] : 0;
-  auto cnt = boost::counting_range<ssize_t>(0, buf.shape[0]);
-  return boost::adaptors::transform(cnt, pairify(buf.ptr, buf.strides[0], stride1));
-  // Be careful that the returned range cannot contain references to dead temporaries.
-}
+#include <pybind11_diagram_utils.h>
 
 double wasserstein_distance(
     Dgm d1, Dgm d2,
diff --git a/src/python/include/pybind11_diagram_utils.h b/src/python/include/pybind11_diagram_utils.h
new file mode 100644
index 00000000..d9627258
--- /dev/null
+++ b/src/python/include/pybind11_diagram_utils.h
@@ -0,0 +1,39 @@
+/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ *    Author(s):       Marc Glisse
+ *
+ *    Copyright (C) 2020 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+
+#include <boost/range/counting_range.hpp>
+#include <boost/range/adaptor/transformed.hpp>
+
+namespace py = pybind11;
+typedef py::array_t<double> Dgm;
+
+// Get m[i,0] and m[i,1] as a pair
+static auto pairify(void* p, ssize_t h, ssize_t w) {
+  return [=](ssize_t i){
+    char* birth = (char*)p + i * h;
+    char* death = birth + w;
+    return std::make_pair(*(double*)birth, *(double*)death);
+  };
+}
+
+inline auto numpy_to_range_of_pairs(py::array_t<double> dgm) {
+  py::buffer_info buf = dgm.request();
+  // shape (n,2) or (0) for empty
+  if((buf.ndim!=2 || buf.shape[1]!=2) && (buf.ndim!=1 || buf.shape[0]!=0))
+    throw std::runtime_error("Diagram must be an array of size n x 2");
+  // In the case of shape (0), avoid reading non-existing strides[1] even if we won't use it.
+  ssize_t stride1 = buf.ndim == 2 ? buf.strides[1] : 0;
+  auto cnt = boost::counting_range<ssize_t>(0, buf.shape[0]);
+  return boost::adaptors::transform(cnt, pairify(buf.ptr, buf.strides[0], stride1));
+  // Be careful that the returned range cannot contain references to dead temporaries.
+}
diff --git a/src/python/setup.py.in b/src/python/setup.py.in
index f968bd59..852da910 100644
--- a/src/python/setup.py.in
+++ b/src/python/setup.py.in
@@ -18,7 +18,8 @@ __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016  Inria"
 __license__ = "MIT"
 
-modules = [@GUDHI_PYTHON_MODULES_TO_COMPILE@]
+cython_modules = [@GUDHI_CYTHON_MODULES@]
+pybind11_modules = [@GUDHI_PYBIND11_MODULES@]
 
 source_dir='@CMAKE_CURRENT_SOURCE_DIR@/gudhi/'
 extra_compile_args=[@GUDHI_PYTHON_EXTRA_COMPILE_ARGS@]
@@ -30,7 +31,7 @@ runtime_library_dirs=[@GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS@]
 
 # Create ext_modules list from module list
 ext_modules = []
-for module in modules:
+for module in cython_modules:
     ext_modules.append(Extension(
         'gudhi.' + module,
         sources = [source_dir + module + '.pyx',],
@@ -55,6 +56,20 @@ ext_modules.append(Extension(
     extra_compile_args=extra_compile_args + [@GUDHI_PYBIND11_EXTRA_COMPILE_ARGS@],
     ))
 
+if "bottleneck" in pybind11_modules:
+    ext_modules.append(Extension(
+        'gudhi.bottleneck',
+        sources = [source_dir + 'bottleneck.cc'],
+        language = 'c++',
+        include_dirs = include_dirs +
+            [pybind11.get_include(False), pybind11.get_include(True)],
+        extra_compile_args=extra_compile_args + [@GUDHI_PYBIND11_EXTRA_COMPILE_ARGS@],
+        extra_link_args=extra_link_args,
+        libraries=libraries,
+        library_dirs=library_dirs,
+        runtime_library_dirs=runtime_library_dirs,
+        ))
+
 setup(
     name = 'gudhi',
     packages=find_packages(), # find_namespace_packages(include=["gudhi*"])