Merge pull request #337 from mglisse/herab

Interface for hera's bottleneck_distance

9 files changed, 97 insertions, 16 deletions

diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 96dd3f6f..1e81cac8 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -130,7 +130,8 @@ if(PYTHONINTERP_FOUND)
     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', ")
+    set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera/wasserstein', ")
+    set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera/bottleneck', ")
     if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
       set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'bottleneck', ")
       set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'nerve_gic', ")
@@ -236,6 +237,7 @@ if(PYTHONINTERP_FOUND)
     file(COPY "gudhi/point_cloud" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
     file(COPY "gudhi/weighted_rips_complex.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
     file(COPY "gudhi/dtm_rips_complex.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
+    file(COPY "gudhi/hera/__init__.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/hera")
 
     add_custom_command(
         OUTPUT gudhi.so
@@ -355,7 +357,9 @@ if(PYTHONINTERP_FOUND)
                COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}"
                ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/bottleneck_basic_example.py")
 
-      add_gudhi_py_test(test_bottleneck_distance)
+      if (PYBIND11_FOUND)
+        add_gudhi_py_test(test_bottleneck_distance)
+      endif()
 
       # Cover complex
       file(COPY ${CMAKE_SOURCE_DIR}/data/points/human.off DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
diff --git a/src/python/doc/bottleneck_distance_user.rst b/src/python/doc/bottleneck_distance_user.rst
index 89da89d3..6c6e08d9 100644
--- a/src/python/doc/bottleneck_distance_user.rst
+++ b/src/python/doc/bottleneck_distance_user.rst
@@ -9,14 +9,23 @@ Definition
 
 .. include:: bottleneck_distance_sum.inc
 
-This implementation is based on ideas from "Geometry Helps in Bottleneck Matching and Related Problems"
-:cite:`DBLP:journals/algorithmica/EfratIK01`. Another relevant publication, although it was not used is
-"Geometry Helps to Compare Persistence Diagrams" :cite:`Kerber:2017:GHC:3047249.3064175`.
+This implementation by François Godi is based on ideas from "Geometry Helps in Bottleneck Matching and Related Problems"
+:cite:`DBLP:journals/algorithmica/EfratIK01` and requires `CGAL <installation.html#cgal>`_ (`GPL v3 </licensing/>`_).
 
-Function
---------
 .. autofunction:: gudhi.bottleneck_distance
 
+This other implementation comes from `Hera
+<https://bitbucket.org/grey_narn/hera/src/master/>`_ (BSD-3-Clause) which is
+based on "Geometry Helps to Compare Persistence Diagrams"
+:cite:`Kerber:2017:GHC:3047249.3064175` by Michael Kerber, Dmitriy
+Morozov, and Arnur Nigmetov.
+
+.. warning::
+   Beware that its approximation allows for a multiplicative error, while the function above uses an additive error.
+
+.. autofunction:: gudhi.hera.bottleneck_distance
+
+
 Distance computation
 --------------------
 
diff --git a/src/python/gudhi/bottleneck.cc b/src/python/gudhi/bottleneck.cc
index 838bf9eb..8a3d669a 100644
--- a/src/python/gudhi/bottleneck.cc
+++ b/src/python/gudhi/bottleneck.cc
@@ -33,7 +33,8 @@ PYBIND11_MODULE(bottleneck, m) {
           py::arg("diagram_1"), py::arg("diagram_2"),
           py::arg("e") = py::none(),
           R"pbdoc(
-    This function returns the point corresponding to a given vertex.
+    Compute the Bottleneck distance between two diagrams.
+    Points at infinity and on the diagonal are supported.
 
     :param diagram_1: The first diagram.
     :type diagram_1: numpy array of shape (m,2)
diff --git a/src/python/gudhi/hera/__init__.py b/src/python/gudhi/hera/__init__.py
new file mode 100644
index 00000000..f70b92b9
--- /dev/null
+++ b/src/python/gudhi/hera/__init__.py
@@ -0,0 +1,7 @@
+from .wasserstein import wasserstein_distance
+from .bottleneck import bottleneck_distance
+
+
+__author__ = "Marc Glisse"
+__copyright__ = "Copyright (C) 2020 Inria"
+__license__ = "MIT"
diff --git a/src/python/gudhi/hera/bottleneck.cc b/src/python/gudhi/hera/bottleneck.cc
new file mode 100644
index 00000000..0cb562ce
--- /dev/null
+++ b/src/python/gudhi/hera/bottleneck.cc
@@ -0,0 +1,54 @@
+/*    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_diagram_utils.h>
+
+#ifdef _MSC_VER
+// https://github.com/grey-narn/hera/issues/3
+// ssize_t is a non-standard type (well, posix)
+using py::ssize_t;
+#endif
+
+#include <bottleneck.h> // Hera
+
+double bottleneck_distance(Dgm d1, Dgm d2, double delta)
+{
+  // 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;
+
+  if (delta == 0)
+    return hera::bottleneckDistExact(diag1, diag2);
+  else
+    return hera::bottleneckDistApprox(diag1, diag2, delta);
+}
+
+PYBIND11_MODULE(bottleneck, m) {
+      m.def("bottleneck_distance", &bottleneck_distance,
+          py::arg("X"), py::arg("Y"),
+          py::arg("delta") = .01,
+          R"pbdoc(
+        Compute the Bottleneck distance between two diagrams.
+        Points at infinity are supported.
+
+        .. note::
+           Points on the diagonal are not supported and must be filtered out before calling this function.
+
+        Parameters:
+            X (n x 2 numpy array): First diagram
+            Y (n x 2 numpy array): Second diagram
+            delta (float): Relative error 1+delta
+
+        Returns:
+            float: (approximate) bottleneck distance d_B(X,Y)
+    )pbdoc");
+}
diff --git a/src/python/gudhi/hera.cc b/src/python/gudhi/hera/wasserstein.cc
index ea80a9a8..1a21f02f 100644
--- a/src/python/gudhi/hera.cc
+++ b/src/python/gudhi/hera/wasserstein.cc
@@ -33,7 +33,7 @@ double wasserstein_distance(
   return hera::wasserstein_dist(diag1, diag2, params);
 }
 
-PYBIND11_MODULE(hera, m) {
+PYBIND11_MODULE(wasserstein, m) {
       m.def("wasserstein_distance", &wasserstein_distance,
           py::arg("X"), py::arg("Y"),
           py::arg("order") = 1,
diff --git a/src/python/include/pybind11_diagram_utils.h b/src/python/include/pybind11_diagram_utils.h
index d9627258..2d5194f4 100644
--- a/src/python/include/pybind11_diagram_utils.h
+++ b/src/python/include/pybind11_diagram_utils.h
@@ -18,8 +18,8 @@ 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){
+static auto pairify(void* p, py::ssize_t h, py::ssize_t w) {
+  return [=](py::ssize_t i){
     char* birth = (char*)p + i * h;
     char* death = birth + w;
     return std::make_pair(*(double*)birth, *(double*)death);
@@ -32,8 +32,8 @@ inline auto numpy_to_range_of_pairs(py::array_t<double> dgm) {
   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]);
+  py::ssize_t stride1 = buf.ndim == 2 ? buf.strides[1] : 0;
+  auto cnt = boost::counting_range<py::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 b9f4e3f0..9b3c7521 100644
--- a/src/python/setup.py.in
+++ b/src/python/setup.py.in
@@ -48,10 +48,12 @@ ext_modules = cythonize(ext_modules)
 
 for module in pybind11_modules:
     my_include_dirs = include_dirs + [pybind11.get_include(False), pybind11.get_include(True)]
-    if module == 'hera':
+    if module == 'hera/wasserstein':
         my_include_dirs = ['@HERA_WASSERSTEIN_INCLUDE_DIR@'] + my_include_dirs
+    elif module == 'hera/bottleneck':
+        my_include_dirs = ['@HERA_BOTTLENECK_INCLUDE_DIR@'] + my_include_dirs
     ext_modules.append(Extension(
-        'gudhi.' + module,
+        'gudhi.' + module.replace('/', '.'),
         sources = [source_dir + module + '.cc'],
         language = 'c++',
         include_dirs = my_include_dirs,
diff --git a/src/python/test/test_bottleneck_distance.py b/src/python/test/test_bottleneck_distance.py
index 70b2abad..6915bea8 100755
--- a/src/python/test/test_bottleneck_distance.py
+++ b/src/python/test/test_bottleneck_distance.py
@@ -9,6 +9,8 @@
 """
 
 import gudhi
+import gudhi.hera
+import pytest
 
 __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016 Inria"
@@ -19,5 +21,7 @@ def test_basic_bottleneck():
     diag1 = [[2.7, 3.7], [9.6, 14.0], [34.2, 34.974], [3.0, float("Inf")]]
     diag2 = [[2.8, 4.45], [9.5, 14.1], [3.2, float("Inf")]]
 
-    assert gudhi.bottleneck_distance(diag1, diag2, 0.1) == 0.8081763781405569
     assert gudhi.bottleneck_distance(diag1, diag2) == 0.75
+    assert gudhi.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, abs=0.1)
+    assert gudhi.hera.bottleneck_distance(diag1, diag2, 0) == 0.75
+    assert gudhi.hera.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, rel=0.1)