diff options
Diffstat (limited to 'src/python/gudhi')
-rw-r--r-- | src/python/gudhi/bottleneck.cc | 51 | ||||
-rw-r--r-- | src/python/gudhi/bottleneck.pyx | 48 | ||||
-rw-r--r-- | src/python/gudhi/hera.cc | 32 |
3 files changed, 52 insertions, 79 deletions
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, |