summaryrefslogtreecommitdiff
path: root/src/python/gudhi
diff options
context:
space:
mode:
Diffstat (limited to 'src/python/gudhi')
-rw-r--r--src/python/gudhi/alpha_complex.pyx49
-rw-r--r--src/python/gudhi/bottleneck.cc15
-rw-r--r--src/python/gudhi/cubical_complex.pyx50
-rw-r--r--src/python/gudhi/dtm_rips_complex.py51
-rw-r--r--src/python/gudhi/hera/__init__.py7
-rw-r--r--src/python/gudhi/hera/bottleneck.cc54
-rw-r--r--src/python/gudhi/hera/wasserstein.cc (renamed from src/python/gudhi/hera.cc)2
-rw-r--r--src/python/gudhi/periodic_cubical_complex.pyx86
-rw-r--r--src/python/gudhi/point_cloud/dtm.py109
-rw-r--r--src/python/gudhi/point_cloud/knn.py4
-rw-r--r--src/python/gudhi/representations/kernel_methods.py41
-rw-r--r--src/python/gudhi/representations/metrics.py56
12 files changed, 407 insertions, 117 deletions
diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx
index d75e374a..a356384d 100644
--- a/src/python/gudhi/alpha_complex.pyx
+++ b/src/python/gudhi/alpha_complex.pyx
@@ -27,11 +27,11 @@ __license__ = "GPL v3"
cdef extern from "Alpha_complex_interface.h" namespace "Gudhi":
cdef cppclass Alpha_complex_interface "Gudhi::alpha_complex::Alpha_complex_interface":
- Alpha_complex_interface(vector[vector[double]] points) nogil except +
+ Alpha_complex_interface(vector[vector[double]] points, bool fast_version) nogil except +
# bool from_file is a workaround for cython to find the correct signature
- Alpha_complex_interface(string off_file, bool from_file) nogil except +
+ Alpha_complex_interface(string off_file, bool fast_version, bool from_file) nogil except +
vector[double] get_point(int vertex) nogil except +
- void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square) nogil except +
+ void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool exact_version, bool default_filtration_value) nogil except +
# AlphaComplex python interface
cdef class AlphaComplex:
@@ -53,10 +53,11 @@ cdef class AlphaComplex:
"""
- cdef Alpha_complex_interface * thisptr
+ cdef Alpha_complex_interface * this_ptr
+ cdef bool exact
# Fake constructor that does nothing but documenting the constructor
- def __init__(self, points=None, off_file=''):
+ def __init__(self, points=None, off_file='', precision='safe'):
"""AlphaComplex constructor.
:param points: A list of points in d-Dimension.
@@ -66,15 +67,21 @@ cdef class AlphaComplex:
:param off_file: An OFF file style name.
:type off_file: string
+
+ :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'.
+ :type precision: string
"""
# The real cython constructor
- def __cinit__(self, points = None, off_file = ''):
+ def __cinit__(self, points = None, off_file = '', precision = 'safe'):
+ assert precision in ['fast', 'safe', 'exact'], "Alpha complex precision can only be 'fast', 'safe' or 'exact'"
+ cdef bool fast = precision == 'fast'
+ self.exact = precision == 'exact'
+
cdef vector[vector[double]] pts
if off_file:
if os.path.isfile(off_file):
- self.thisptr = new Alpha_complex_interface(
- off_file.encode('utf-8'), True)
+ self.this_ptr = new Alpha_complex_interface(off_file.encode('utf-8'), fast, True)
else:
print("file " + off_file + " not found.")
else:
@@ -83,17 +90,16 @@ cdef class AlphaComplex:
points=[]
pts = points
with nogil:
- self.thisptr = new Alpha_complex_interface(pts)
-
+ self.this_ptr = new Alpha_complex_interface(pts, fast)
def __dealloc__(self):
- if self.thisptr != NULL:
- del self.thisptr
+ if self.this_ptr != NULL:
+ del self.this_ptr
def __is_defined(self):
"""Returns true if AlphaComplex pointer is not NULL.
"""
- return self.thisptr != NULL
+ return self.this_ptr != NULL
def get_point(self, vertex):
"""This function returns the point corresponding to a given vertex.
@@ -103,21 +109,24 @@ cdef class AlphaComplex:
:rtype: list of float
:returns: the point.
"""
- return self.thisptr.get_point(vertex)
+ return self.this_ptr.get_point(vertex)
- def create_simplex_tree(self, max_alpha_square = float('inf')):
+ def create_simplex_tree(self, max_alpha_square = float('inf'), default_filtration_value = False):
"""
- :param max_alpha_square: The maximum alpha square threshold the
- simplices shall not exceed. Default is set to infinity, and
- there is very little point using anything else since it does
- not save time.
+ :param max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to
+ infinity, and there is very little point using anything else since it does not save time.
:type max_alpha_square: float
+ :param default_filtration_value: Set this value to `True` if filtration values are not needed to be computed
+ (will be set to `NaN`). Default value is `False` (which means compute the filtration values).
+ :type default_filtration_value: bool
:returns: A simplex tree created from the Delaunay Triangulation.
:rtype: SimplexTree
"""
stree = SimplexTree()
cdef double mas = max_alpha_square
cdef intptr_t stree_int_ptr=stree.thisptr
+ cdef bool compute_filtration = default_filtration_value == True
with nogil:
- self.thisptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr, mas)
+ self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr,
+ mas, self.exact, compute_filtration)
return stree
diff --git a/src/python/gudhi/bottleneck.cc b/src/python/gudhi/bottleneck.cc
index 732cb9a8..8a3d669a 100644
--- a/src/python/gudhi/bottleneck.cc
+++ b/src/python/gudhi/bottleneck.cc
@@ -12,24 +12,29 @@
#include <pybind11_diagram_utils.h>
-double bottleneck(Dgm d1, Dgm d2, double epsilon)
+// For compatibility with older versions, we want to support e=None.
+// In C++17, the recommended way is std::optional<double>.
+double bottleneck(Dgm d1, Dgm d2, py::object epsilon)
{
+ double e = (std::numeric_limits<double>::min)();
+ if (!epsilon.is_none()) e = epsilon.cast<double>();
// 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);
+ return Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, e);
}
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)(),
+ 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)
@@ -42,7 +47,7 @@ PYBIND11_MODULE(bottleneck, m) {
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.
+ Thus, by default (`e=None`), `e` is the smallest positive double.
:type e: float
:rtype: float
:returns: the bottleneck distance.
diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx
index ca979eda..28fbe3af 100644
--- a/src/python/gudhi/cubical_complex.pyx
+++ b/src/python/gudhi/cubical_complex.pyx
@@ -27,20 +27,20 @@ __license__ = "MIT"
cdef extern from "Cubical_complex_interface.h" namespace "Gudhi":
cdef cppclass Bitmap_cubical_complex_base_interface "Gudhi::Cubical_complex::Cubical_complex_interface<>":
- Bitmap_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells)
- Bitmap_cubical_complex_base_interface(string perseus_file)
- int num_simplices()
- int dimension()
+ Bitmap_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells) nogil
+ Bitmap_cubical_complex_base_interface(string perseus_file) nogil
+ int num_simplices() nogil
+ int dimension() nogil
cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
cdef cppclass Cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Cubical_complex::Cubical_complex_interface<>>":
- Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max)
- void compute_persistence(int homology_coeff_field, double min_persistence)
- vector[pair[int, pair[double, double]]] get_persistence()
- vector[vector[int]] cofaces_of_cubical_persistence_pairs()
- vector[int] betti_numbers()
- vector[int] persistent_betti_numbers(double from_value, double to_value)
- vector[pair[double,double]] intervals_in_dimension(int dimension)
+ Cubical_complex_persistence_interface(Bitmap_cubical_complex_base_interface * st, bool persistence_dim_max) nogil
+ void compute_persistence(int homology_coeff_field, double min_persistence) nogil
+ vector[pair[int, pair[double, double]]] get_persistence() nogil
+ vector[vector[int]] cofaces_of_cubical_persistence_pairs() nogil
+ vector[int] betti_numbers() nogil
+ vector[int] persistent_betti_numbers(double from_value, double to_value) nogil
+ vector[pair[double,double]] intervals_in_dimension(int dimension) nogil
# CubicalComplex python interface
cdef class CubicalComplex:
@@ -80,7 +80,7 @@ cdef class CubicalComplex:
perseus_file=''):
if ((dimensions is not None) and (top_dimensional_cells is not None)
and (perseus_file == '')):
- self.thisptr = new Bitmap_cubical_complex_base_interface(dimensions, top_dimensional_cells)
+ self._construct_from_cells(dimensions, top_dimensional_cells)
elif ((dimensions is None) and (top_dimensional_cells is not None)
and (perseus_file == '')):
top_dimensional_cells = np.array(top_dimensional_cells,
@@ -88,11 +88,11 @@ cdef class CubicalComplex:
order = 'F')
dimensions = top_dimensional_cells.shape
top_dimensional_cells = top_dimensional_cells.ravel(order='F')
- self.thisptr = new Bitmap_cubical_complex_base_interface(dimensions, top_dimensional_cells)
+ self._construct_from_cells(dimensions, top_dimensional_cells)
elif ((dimensions is None) and (top_dimensional_cells is None)
and (perseus_file != '')):
if os.path.isfile(perseus_file):
- self.thisptr = new Bitmap_cubical_complex_base_interface(perseus_file.encode('utf-8'))
+ self._construct_from_file(perseus_file.encode('utf-8'))
else:
raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT),
perseus_file)
@@ -101,6 +101,14 @@ cdef class CubicalComplex:
"top_dimensional_cells or from a Perseus-style file name.",
file=sys.stderr)
+ def _construct_from_cells(self, vector[unsigned] dimensions, vector[double] top_dimensional_cells):
+ with nogil:
+ self.thisptr = new Bitmap_cubical_complex_base_interface(dimensions, top_dimensional_cells)
+
+ def _construct_from_file(self, string filename):
+ with nogil:
+ self.thisptr = new Bitmap_cubical_complex_base_interface(filename)
+
def __dealloc__(self):
if self.thisptr != NULL:
del self.thisptr
@@ -151,8 +159,11 @@ cdef class CubicalComplex:
if self.pcohptr != NULL:
del self.pcohptr
assert self.__is_defined()
- self.pcohptr = new Cubical_complex_persistence_interface(self.thisptr, True)
- self.pcohptr.compute_persistence(homology_coeff_field, min_persistence)
+ cdef int field = homology_coeff_field
+ cdef double minp = min_persistence
+ with nogil:
+ self.pcohptr = new Cubical_complex_persistence_interface(self.thisptr, 1)
+ self.pcohptr.compute_persistence(field, minp)
def persistence(self, homology_coeff_field=11, min_persistence=0):
"""This function computes and returns the persistence of the complex.
@@ -204,19 +215,20 @@ cdef class CubicalComplex:
cdef vector[vector[int]] persistence_result
output = [[],[]]
- persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs()
+ with nogil:
+ persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs()
pr = np.array(persistence_result)
ess_ind = np.argwhere(pr[:,2] == -1)[:,0]
ess = pr[ess_ind]
- max_h = max(ess[:,0])+1
+ max_h = max(ess[:,0])+1 if len(ess) > 0 else 0
for h in range(max_h):
hidxs = np.argwhere(ess[:,0] == h)[:,0]
output[1].append(ess[hidxs][:,1])
reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind)
reg = pr[reg_ind]
- max_h = max(reg[:,0])+1
+ max_h = max(reg[:,0])+1 if len(reg) > 0 else 0
for h in range(max_h):
hidxs = np.argwhere(reg[:,0] == h)[:,0]
output[0].append(reg[hidxs][:,1:])
diff --git a/src/python/gudhi/dtm_rips_complex.py b/src/python/gudhi/dtm_rips_complex.py
new file mode 100644
index 00000000..63c9b138
--- /dev/null
+++ b/src/python/gudhi/dtm_rips_complex.py
@@ -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): Yuichi Ike, Raphaël Tinarrage
+#
+# Copyright (C) 2020 Inria, Copyright (C) 2020 FUjitsu Laboratories Ltd.
+#
+# Modification(s):
+# - YYYY/MM Author: Description of the modification
+
+
+from gudhi.weighted_rips_complex import WeightedRipsComplex
+from gudhi.point_cloud.dtm import DistanceToMeasure
+from scipy.spatial.distance import cdist
+
+class DTMRipsComplex(WeightedRipsComplex):
+ """
+ Class to generate a DTM Rips complex from a distance matrix or a point set,
+ in the way described in :cite:`dtmfiltrations`.
+ Remark that all the filtration values are doubled compared to the definition in the paper
+ for the consistency with RipsComplex.
+ :Requires: `SciPy <installation.html#scipy>`_
+ """
+ def __init__(self,
+ points=None,
+ distance_matrix=None,
+ k=1,
+ q=2,
+ max_filtration=float('inf')):
+ """
+ Args:
+ points (numpy.ndarray): array of points.
+ distance_matrix (numpy.ndarray): full distance matrix.
+ k (int): number of neighbors for the computation of DTM. Defaults to 1, which is equivalent to the usual Rips complex.
+ q (float): order used to compute the distance to measure. Defaults to 2.
+ max_filtration (float): specifies the maximal filtration value to be considered.
+ """
+ if distance_matrix is None:
+ if points is None:
+ # Empty Rips construction
+ points=[]
+ distance_matrix = cdist(points,points)
+ self.distance_matrix = distance_matrix
+
+ # TODO: address the error when k is too large
+ if k <= 1:
+ self.weights = [0] * len(distance_matrix)
+ else:
+ dtm = DistanceToMeasure(k, q=q, metric="precomputed")
+ self.weights = dtm.fit_transform(distance_matrix)
+ self.max_filtration = max_filtration
+
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/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx
index 06309772..d353d2af 100644
--- a/src/python/gudhi/periodic_cubical_complex.pyx
+++ b/src/python/gudhi/periodic_cubical_complex.pyx
@@ -24,20 +24,20 @@ __license__ = "MIT"
cdef extern from "Cubical_complex_interface.h" namespace "Gudhi":
cdef cppclass Periodic_cubical_complex_base_interface "Gudhi::Cubical_complex::Cubical_complex_interface<Gudhi::cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double>>":
- Periodic_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells, vector[bool] periodic_dimensions)
- Periodic_cubical_complex_base_interface(string perseus_file)
- int num_simplices()
- int dimension()
+ Periodic_cubical_complex_base_interface(vector[unsigned] dimensions, vector[double] top_dimensional_cells, vector[bool] periodic_dimensions) nogil
+ Periodic_cubical_complex_base_interface(string perseus_file) nogil
+ int num_simplices() nogil
+ int dimension() nogil
cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
cdef cppclass Periodic_cubical_complex_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Cubical_complex::Cubical_complex_interface<Gudhi::cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double>>>":
- Periodic_cubical_complex_persistence_interface(Periodic_cubical_complex_base_interface * st, bool persistence_dim_max)
- void compute_persistence(int homology_coeff_field, double min_persistence)
- vector[pair[int, pair[double, double]]] get_persistence()
- vector[vector[int]] cofaces_of_cubical_persistence_pairs()
- vector[int] betti_numbers()
- vector[int] persistent_betti_numbers(double from_value, double to_value)
- vector[pair[double,double]] intervals_in_dimension(int dimension)
+ Periodic_cubical_complex_persistence_interface(Periodic_cubical_complex_base_interface * st, bool persistence_dim_max) nogil
+ void compute_persistence(int homology_coeff_field, double min_persistence) nogil
+ vector[pair[int, pair[double, double]]] get_persistence() nogil
+ vector[vector[int]] cofaces_of_cubical_persistence_pairs() nogil
+ vector[int] betti_numbers() nogil
+ vector[int] persistent_betti_numbers(double from_value, double to_value) nogil
+ vector[pair[double,double]] intervals_in_dimension(int dimension) nogil
# PeriodicCubicalComplex python interface
cdef class PeriodicCubicalComplex:
@@ -81,9 +81,7 @@ cdef class PeriodicCubicalComplex:
periodic_dimensions=None, perseus_file=''):
if ((dimensions is not None) and (top_dimensional_cells is not None)
and (periodic_dimensions is not None) and (perseus_file == '')):
- self.thisptr = new Periodic_cubical_complex_base_interface(dimensions,
- top_dimensional_cells,
- periodic_dimensions)
+ self._construct_from_cells(dimensions, top_dimensional_cells, periodic_dimensions)
elif ((dimensions is None) and (top_dimensional_cells is not None)
and (periodic_dimensions is not None) and (perseus_file == '')):
top_dimensional_cells = np.array(top_dimensional_cells,
@@ -91,13 +89,11 @@ cdef class PeriodicCubicalComplex:
order = 'F')
dimensions = top_dimensional_cells.shape
top_dimensional_cells = top_dimensional_cells.ravel(order='F')
- self.thisptr = new Periodic_cubical_complex_base_interface(dimensions,
- top_dimensional_cells,
- periodic_dimensions)
+ self._construct_from_cells(dimensions, top_dimensional_cells, periodic_dimensions)
elif ((dimensions is None) and (top_dimensional_cells is None)
and (periodic_dimensions is None) and (perseus_file != '')):
if os.path.isfile(perseus_file):
- self.thisptr = new Periodic_cubical_complex_base_interface(perseus_file.encode('utf-8'))
+ self._construct_from_file(perseus_file.encode('utf-8'))
else:
print("file " + perseus_file + " not found.", file=sys.stderr)
else:
@@ -106,6 +102,14 @@ cdef class PeriodicCubicalComplex:
"top_dimensional_cells and periodic_dimensions or from "
"a Perseus-style file name.", file=sys.stderr)
+ def _construct_from_cells(self, vector[unsigned] dimensions, vector[double] top_dimensional_cells, vector[bool] periodic_dimensions):
+ with nogil:
+ self.thisptr = new Periodic_cubical_complex_base_interface(dimensions, top_dimensional_cells, periodic_dimensions)
+
+ def _construct_from_file(self, string filename):
+ with nogil:
+ self.thisptr = new Periodic_cubical_complex_base_interface(filename)
+
def __dealloc__(self):
if self.thisptr != NULL:
del self.thisptr
@@ -156,8 +160,11 @@ cdef class PeriodicCubicalComplex:
if self.pcohptr != NULL:
del self.pcohptr
assert self.__is_defined()
- self.pcohptr = new Periodic_cubical_complex_persistence_interface(self.thisptr, True)
- self.pcohptr.compute_persistence(homology_coeff_field, min_persistence)
+ cdef int field = homology_coeff_field
+ cdef double minp = min_persistence
+ with nogil:
+ self.pcohptr = new Periodic_cubical_complex_persistence_interface(self.thisptr, 1)
+ self.pcohptr.compute_persistence(field, minp)
def persistence(self, homology_coeff_field=11, min_persistence=0):
"""This function computes and returns the persistence of the complex.
@@ -204,28 +211,27 @@ cdef class PeriodicCubicalComplex:
The indices of the arrays in the list correspond to the homological dimensions, and the
integers of each row in each array correspond to: (index of positive top-dimensional cell).
"""
+ assert self.pcohptr != NULL, "compute_persistence() must be called before cofaces_of_persistence_pairs()"
cdef vector[vector[int]] persistence_result
- if self.pcohptr != NULL:
- output = [[],[]]
+
+ output = [[],[]]
+ with nogil:
persistence_result = self.pcohptr.cofaces_of_cubical_persistence_pairs()
- pr = np.array(persistence_result)
-
- ess_ind = np.argwhere(pr[:,2] == -1)[:,0]
- ess = pr[ess_ind]
- max_h = max(ess[:,0])+1
- for h in range(max_h):
- hidxs = np.argwhere(ess[:,0] == h)[:,0]
- output[1].append(ess[hidxs][:,1])
-
- reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind)
- reg = pr[reg_ind]
- max_h = max(reg[:,0])+1
- for h in range(max_h):
- hidxs = np.argwhere(reg[:,0] == h)[:,0]
- output[0].append(reg[hidxs][:,1:])
- else:
- print("cofaces_of_persistence_pairs function requires persistence function"
- " to be launched first.")
+ pr = np.array(persistence_result)
+
+ ess_ind = np.argwhere(pr[:,2] == -1)[:,0]
+ ess = pr[ess_ind]
+ max_h = max(ess[:,0])+1 if len(ess) > 0 else 0
+ for h in range(max_h):
+ hidxs = np.argwhere(ess[:,0] == h)[:,0]
+ output[1].append(ess[hidxs][:,1])
+
+ reg_ind = np.setdiff1d(np.array(range(len(pr))), ess_ind)
+ reg = pr[reg_ind]
+ max_h = max(reg[:,0])+1 if len(reg) > 0 else 0
+ for h in range(max_h):
+ hidxs = np.argwhere(reg[:,0] == h)[:,0]
+ output[0].append(reg[hidxs][:,1:])
return output
def betti_numbers(self):
diff --git a/src/python/gudhi/point_cloud/dtm.py b/src/python/gudhi/point_cloud/dtm.py
index 13e16d24..55ac58e6 100644
--- a/src/python/gudhi/point_cloud/dtm.py
+++ b/src/python/gudhi/point_cloud/dtm.py
@@ -8,6 +8,7 @@
# - YYYY/MM Author: Description of the modification
from .knn import KNearestNeighbors
+import numpy as np
__author__ = "Marc Glisse"
__copyright__ = "Copyright (C) 2020 Inria"
@@ -68,3 +69,111 @@ class DistanceToMeasure:
# We compute too many powers, 1/p in knn then q in dtm, 1/q in dtm then q or some log in the caller.
# Add option to skip the final root?
return dtm
+
+
+class DTMDensity:
+ """
+ Density estimator based on the distance to the empirical measure defined by a point set, as defined
+ in :cite:`dtmdensity`. Note that this implementation only renormalizes when asked, and the renormalization
+ only works for a Euclidean metric, so in other cases the total measure may not be 1.
+
+ .. note:: When the dimension is high, using it as an exponent can quickly lead to under- or overflows.
+ We recommend using a small fixed value instead in those cases, even if it won't have the same nice
+ theoretical properties as the dimension.
+ """
+
+ def __init__(self, k=None, weights=None, q=None, dim=None, normalize=False, n_samples=None, **kwargs):
+ """
+ Args:
+ k (int): number of neighbors (possibly including the point itself). Optional if it can be guessed
+ from weights or metric="neighbors".
+ weights (numpy.array): weights of each of the k neighbors, optional. They are supposed to sum to 1.
+ q (float): order used to compute the distance to measure. Defaults to dim.
+ dim (float): final exponent representing the dimension. Defaults to the dimension, and must be specified
+ when the dimension cannot be read from the input (metric is "neighbors" or "precomputed").
+ normalize (bool): normalize the density so it corresponds to a probability measure on ℝᵈ.
+ Only available for the Euclidean metric, defaults to False.
+ n_samples (int): number of sample points used for fitting. Only needed if `normalize` is True and
+ metric is "neighbors".
+ kwargs: same parameters as :class:`~gudhi.point_cloud.knn.KNearestNeighbors`, except that
+ metric="neighbors" means that :func:`transform` expects an array with the distances to
+ the k nearest neighbors.
+ """
+ if weights is None:
+ self.k = k
+ if k is None:
+ assert kwargs.get("metric") == "neighbors", 'Must specify k or weights, unless metric is "neighbors"'
+ self.weights = None
+ else:
+ self.weights = np.full(k, 1.0 / k)
+ else:
+ self.weights = weights
+ self.k = len(weights)
+ assert k is None or k == self.k, "k differs from the length of weights"
+ self.q = q
+ self.dim = dim
+ self.params = kwargs
+ self.normalize = normalize
+ self.n_samples = n_samples
+
+ def fit_transform(self, X, y=None):
+ return self.fit(X).transform(X)
+
+ def fit(self, X, y=None):
+ """
+ Args:
+ X (numpy.array): coordinates for mass points.
+ """
+ if self.params.setdefault("metric", "euclidean") != "neighbors":
+ self.knn = KNearestNeighbors(
+ self.k, return_index=False, return_distance=True, sort_results=False, **self.params
+ )
+ self.knn.fit(X)
+ if self.params["metric"] != "precomputed":
+ self.n_samples = len(X)
+ return self
+
+ def transform(self, X):
+ """
+ Args:
+ X (numpy.array): coordinates for query points, or distance matrix if metric is "precomputed",
+ or distances to the k nearest neighbors if metric is "neighbors" (if the array has more
+ than k columns, the remaining ones are ignored).
+ """
+ q = self.q
+ dim = self.dim
+ if dim is None:
+ assert self.params["metric"] not in {
+ "neighbors",
+ "precomputed",
+ }, "dim not specified and cannot guess the dimension"
+ dim = len(X[0])
+ if q is None:
+ q = dim
+ k = self.k
+ weights = self.weights
+ if self.params["metric"] == "neighbors":
+ distances = np.asarray(X)
+ if weights is None:
+ k = distances.shape[1]
+ weights = np.full(k, 1.0 / k)
+ else:
+ distances = distances[:, :k]
+ else:
+ distances = self.knn.transform(X)
+ distances = distances ** q
+ dtm = (distances * weights).sum(-1)
+ if self.normalize:
+ dtm /= (np.arange(1, k + 1) ** (q / dim) * weights).sum()
+ density = dtm ** (-dim / q)
+ if self.normalize:
+ import math
+
+ if self.params["metric"] == "precomputed":
+ self.n_samples = len(X[0])
+ # Volume of d-ball
+ Vd = math.pi ** (dim / 2) / math.gamma(dim / 2 + 1)
+ density /= self.n_samples * Vd
+ return density
+ # We compute too many powers, 1/p in knn then q in dtm, d/q in dtm then whatever in the caller.
+ # Add option to skip the final root?
diff --git a/src/python/gudhi/point_cloud/knn.py b/src/python/gudhi/point_cloud/knn.py
index 86008bc3..4652fe80 100644
--- a/src/python/gudhi/point_cloud/knn.py
+++ b/src/python/gudhi/point_cloud/knn.py
@@ -306,6 +306,10 @@ class KNearestNeighbors:
if self.params["implementation"] == "ckdtree":
qargs = {key: val for key, val in self.params.items() if key in {"p", "eps", "n_jobs"}}
distances, neighbors = self.kdtree.query(X, k=self.k, **qargs)
+ if k == 1:
+ # SciPy decided to squeeze the last dimension for k=1
+ distances = distances[:, None]
+ neighbors = neighbors[:, None]
if self.return_index:
if self.return_distance:
return neighbors, distances
diff --git a/src/python/gudhi/representations/kernel_methods.py b/src/python/gudhi/representations/kernel_methods.py
index 596f4f07..23fd23c7 100644
--- a/src/python/gudhi/representations/kernel_methods.py
+++ b/src/python/gudhi/representations/kernel_methods.py
@@ -10,7 +10,7 @@
import numpy as np
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.metrics import pairwise_distances, pairwise_kernels
-from .metrics import SlicedWassersteinDistance, PersistenceFisherDistance, _sklearn_wrapper, pairwise_persistence_diagram_distances, _sliced_wasserstein_distance, _persistence_fisher_distance
+from .metrics import SlicedWassersteinDistance, PersistenceFisherDistance, _sklearn_wrapper, _pairwise, pairwise_persistence_diagram_distances, _sliced_wasserstein_distance, _persistence_fisher_distance
from .preprocessing import Padding
#############################################
@@ -60,7 +60,7 @@ def _persistence_scale_space_kernel(D1, D2, kernel_approx=None, bandwidth=1.):
weight_pss = lambda x: 1 if x[1] >= x[0] else -1
return 0.5 * _persistence_weighted_gaussian_kernel(DD1, DD2, weight=weight_pss, kernel_approx=kernel_approx, bandwidth=bandwidth)
-def pairwise_persistence_diagram_kernels(X, Y=None, kernel="sliced_wasserstein", **kwargs):
+def pairwise_persistence_diagram_kernels(X, Y=None, kernel="sliced_wasserstein", n_jobs=None, **kwargs):
"""
This function computes the kernel matrix between two lists of persistence diagrams given as numpy arrays of shape (nx2).
@@ -68,38 +68,41 @@ def pairwise_persistence_diagram_kernels(X, Y=None, kernel="sliced_wasserstein",
X (list of n numpy arrays of shape (numx2)): first list of persistence diagrams.
Y (list of m numpy arrays of shape (numx2)): second list of persistence diagrams (optional). If None, pairwise kernel values are computed from the first list only.
kernel: kernel to use. It can be either a string ("sliced_wasserstein", "persistence_scale_space", "persistence_weighted_gaussian", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs. If it is a function, make sure that it is symmetric.
+ n_jobs (int): number of jobs to use for the computation. This uses joblib.Parallel(prefer="threads"), so kernels that do not release the GIL may not scale unless run inside a `joblib.parallel_backend <https://joblib.readthedocs.io/en/latest/parallel.html#joblib.parallel_backend>`_ block.
**kwargs: optional keyword parameters. Any further parameters are passed directly to the kernel function. See the docs of the various kernel classes in this module.
Returns:
numpy array of shape (nxm): kernel matrix.
"""
XX = np.reshape(np.arange(len(X)), [-1,1])
- YY = None if Y is None else np.reshape(np.arange(len(Y)), [-1,1])
+ YY = None if Y is None or Y is X else np.reshape(np.arange(len(Y)), [-1,1])
if kernel == "sliced_wasserstein":
- return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="sliced_wasserstein", num_directions=kwargs["num_directions"]) / kwargs["bandwidth"])
+ return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="sliced_wasserstein", num_directions=kwargs["num_directions"], n_jobs=n_jobs) / kwargs["bandwidth"])
elif kernel == "persistence_fisher":
- return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="persistence_fisher", kernel_approx=kwargs["kernel_approx"], bandwidth=kwargs["bandwidth"]) / kwargs["bandwidth_fisher"])
+ return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="persistence_fisher", kernel_approx=kwargs["kernel_approx"], bandwidth=kwargs["bandwidth"], n_jobs=n_jobs) / kwargs["bandwidth_fisher"])
elif kernel == "persistence_scale_space":
- return pairwise_kernels(XX, YY, metric=_sklearn_wrapper(_persistence_scale_space_kernel, X, Y, **kwargs))
+ return _pairwise(pairwise_kernels, False, XX, YY, metric=_sklearn_wrapper(_persistence_scale_space_kernel, X, Y, **kwargs), n_jobs=n_jobs)
elif kernel == "persistence_weighted_gaussian":
- return pairwise_kernels(XX, YY, metric=_sklearn_wrapper(_persistence_weighted_gaussian_kernel, X, Y, **kwargs))
+ return _pairwise(pairwise_kernels, False, XX, YY, metric=_sklearn_wrapper(_persistence_weighted_gaussian_kernel, X, Y, **kwargs), n_jobs=n_jobs)
else:
- return pairwise_kernels(XX, YY, metric=_sklearn_wrapper(metric, **kwargs))
+ return _pairwise(pairwise_kernels, False, XX, YY, metric=_sklearn_wrapper(metric, **kwargs), n_jobs=n_jobs)
class SlicedWassersteinKernel(BaseEstimator, TransformerMixin):
"""
This is a class for computing the sliced Wasserstein kernel matrix from a list of persistence diagrams. The sliced Wasserstein kernel is computed by exponentiating the corresponding sliced Wasserstein distance with a Gaussian kernel. See http://proceedings.mlr.press/v70/carriere17a.html for more details.
"""
- def __init__(self, num_directions=10, bandwidth=1.0):
+ def __init__(self, num_directions=10, bandwidth=1.0, n_jobs=None):
"""
Constructor for the SlicedWassersteinKernel class.
Parameters:
bandwidth (double): bandwidth of the Gaussian kernel applied to the sliced Wasserstein distance (default 1.).
num_directions (int): number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the kernel computation (default 10).
+ n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_kernels` for details.
"""
self.bandwidth = bandwidth
self.num_directions = num_directions
+ self.n_jobs = n_jobs
def fit(self, X, y=None):
"""
@@ -122,7 +125,7 @@ class SlicedWassersteinKernel(BaseEstimator, TransformerMixin):
Returns:
numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise sliced Wasserstein kernel values.
"""
- return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="sliced_wasserstein", bandwidth=self.bandwidth, num_directions=self.num_directions)
+ return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="sliced_wasserstein", bandwidth=self.bandwidth, num_directions=self.num_directions, n_jobs=self.n_jobs)
def __call__(self, diag1, diag2):
"""
@@ -141,7 +144,7 @@ class PersistenceWeightedGaussianKernel(BaseEstimator, TransformerMixin):
"""
This is a class for computing the persistence weighted Gaussian kernel matrix from a list of persistence diagrams. The persistence weighted Gaussian kernel is computed by convolving the persistence diagram points with weighted Gaussian kernels. See http://proceedings.mlr.press/v48/kusano16.html for more details.
"""
- def __init__(self, bandwidth=1., weight=lambda x: 1, kernel_approx=None):
+ def __init__(self, bandwidth=1., weight=lambda x: 1, kernel_approx=None, n_jobs=None):
"""
Constructor for the PersistenceWeightedGaussianKernel class.
@@ -149,9 +152,11 @@ class PersistenceWeightedGaussianKernel(BaseEstimator, TransformerMixin):
bandwidth (double): bandwidth of the Gaussian kernel with which persistence diagrams will be convolved (default 1.)
weight (function): weight function for the persistence diagram points (default constant function, ie lambda x: 1). This function must be defined on 2D points, ie lists or numpy arrays of the form [p_x,p_y].
kernel_approx (class): kernel approximation class used to speed up computation (default None). Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
+ n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_kernels` for details.
"""
self.bandwidth, self.weight = bandwidth, weight
self.kernel_approx = kernel_approx
+ self.n_jobs = n_jobs
def fit(self, X, y=None):
"""
@@ -174,7 +179,7 @@ class PersistenceWeightedGaussianKernel(BaseEstimator, TransformerMixin):
Returns:
numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence weighted Gaussian kernel values.
"""
- return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_weighted_gaussian", bandwidth=self.bandwidth, weight=self.weight, kernel_approx=self.kernel_approx)
+ return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_weighted_gaussian", bandwidth=self.bandwidth, weight=self.weight, kernel_approx=self.kernel_approx, n_jobs=self.n_jobs)
def __call__(self, diag1, diag2):
"""
@@ -193,15 +198,17 @@ class PersistenceScaleSpaceKernel(BaseEstimator, TransformerMixin):
"""
This is a class for computing the persistence scale space kernel matrix from a list of persistence diagrams. The persistence scale space kernel is computed by adding the symmetric to the diagonal of each point in each persistence diagram, with negative weight, and then convolving the points with a Gaussian kernel. See https://www.cv-foundation.org/openaccess/content_cvpr_2015/papers/Reininghaus_A_Stable_Multi-Scale_2015_CVPR_paper.pdf for more details.
"""
- def __init__(self, bandwidth=1., kernel_approx=None):
+ def __init__(self, bandwidth=1., kernel_approx=None, n_jobs=None):
"""
Constructor for the PersistenceScaleSpaceKernel class.
Parameters:
bandwidth (double): bandwidth of the Gaussian kernel with which persistence diagrams will be convolved (default 1.)
kernel_approx (class): kernel approximation class used to speed up computation (default None). Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
+ n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_kernels` for details.
"""
self.bandwidth, self.kernel_approx = bandwidth, kernel_approx
+ self.n_jobs = n_jobs
def fit(self, X, y=None):
"""
@@ -224,7 +231,7 @@ class PersistenceScaleSpaceKernel(BaseEstimator, TransformerMixin):
Returns:
numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence scale space kernel values.
"""
- return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_scale_space", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)
+ return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_scale_space", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx, n_jobs=self.n_jobs)
def __call__(self, diag1, diag2):
"""
@@ -243,7 +250,7 @@ class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
"""
This is a class for computing the persistence Fisher kernel matrix from a list of persistence diagrams. The persistence Fisher kernel is computed by exponentiating the corresponding persistence Fisher distance with a Gaussian kernel. See papers.nips.cc/paper/8205-persistence-fisher-kernel-a-riemannian-manifold-kernel-for-persistence-diagrams for more details.
"""
- def __init__(self, bandwidth_fisher=1., bandwidth=1., kernel_approx=None):
+ def __init__(self, bandwidth_fisher=1., bandwidth=1., kernel_approx=None, n_jobs=None):
"""
Constructor for the PersistenceFisherKernel class.
@@ -251,9 +258,11 @@ class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
bandwidth (double): bandwidth of the Gaussian kernel applied to the persistence Fisher distance (default 1.).
bandwidth_fisher (double): bandwidth of the Gaussian kernel used to turn persistence diagrams into probability distributions by PersistenceFisherDistance class (default 1.).
kernel_approx (class): kernel approximation class used to speed up computation (default None). Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
+ n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_kernels` for details.
"""
self.bandwidth = bandwidth
self.bandwidth_fisher, self.kernel_approx = bandwidth_fisher, kernel_approx
+ self.n_jobs = n_jobs
def fit(self, X, y=None):
"""
@@ -276,7 +285,7 @@ class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
Returns:
numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence Fisher kernel values.
"""
- return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_fisher", bandwidth=self.bandwidth, bandwidth_fisher=self.bandwidth_fisher, kernel_approx=self.kernel_approx)
+ return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_fisher", bandwidth=self.bandwidth, bandwidth_fisher=self.bandwidth_fisher, kernel_approx=self.kernel_approx, n_jobs=self.n_jobs)
def __call__(self, diag1, diag2):
"""
diff --git a/src/python/gudhi/representations/metrics.py b/src/python/gudhi/representations/metrics.py
index 8a32f7e9..cf2e0879 100644
--- a/src/python/gudhi/representations/metrics.py
+++ b/src/python/gudhi/representations/metrics.py
@@ -12,6 +12,7 @@ from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.metrics import pairwise_distances
from gudhi.hera import wasserstein_distance as hera_wasserstein_distance
from .preprocessing import Padding
+from joblib import Parallel, delayed
#############################################
# Metrics ###################################
@@ -116,6 +117,20 @@ def _persistence_fisher_distance(D1, D2, kernel_approx=None, bandwidth=1.):
vectorj = vectorj/vectorj_sum
return np.arccos( min(np.dot(np.sqrt(vectori), np.sqrt(vectorj)), 1.) )
+def _pairwise(fallback, skipdiag, X, Y, metric, n_jobs):
+ if Y is not None:
+ return fallback(X, Y, metric=metric, n_jobs=n_jobs)
+ triu = np.triu_indices(len(X), k=skipdiag)
+ tril = (triu[1], triu[0])
+ par = Parallel(n_jobs=n_jobs, prefer="threads")
+ d = par(delayed(metric)([triu[0][i]], [triu[1][i]]) for i in range(len(triu[0])))
+ m = np.empty((len(X), len(X)))
+ m[triu] = d
+ m[tril] = d
+ if skipdiag:
+ np.fill_diagonal(m, 0)
+ return m
+
def _sklearn_wrapper(metric, X, Y, **kwargs):
"""
This function is a wrapper for any metric between two persistence diagrams that takes two numpy arrays of shapes (nx2) and (mx2) as arguments.
@@ -134,7 +149,7 @@ PAIRWISE_DISTANCE_FUNCTIONS = {
"persistence_fisher": _persistence_fisher_distance,
}
-def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwargs):
+def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", n_jobs=None, **kwargs):
"""
This function computes the distance matrix between two lists of persistence diagrams given as numpy arrays of shape (nx2).
@@ -142,48 +157,51 @@ def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwa
X (list of n numpy arrays of shape (numx2)): first list of persistence diagrams.
Y (list of m numpy arrays of shape (numx2)): second list of persistence diagrams (optional). If None, pairwise distances are computed from the first list only.
metric: distance to use. It can be either a string ("sliced_wasserstein", "wasserstein", "hera_wasserstein" (Wasserstein distance computed with Hera---note that Hera is also used for the default option "wasserstein"), "pot_wasserstein" (Wasserstein distance computed with POT), "bottleneck", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs. If it is a function, make sure that it is symmetric and that it outputs 0 if called on the same two arrays.
+ n_jobs (int): number of jobs to use for the computation. This uses joblib.Parallel(prefer="threads"), so metrics that do not release the GIL may not scale unless run inside a `joblib.parallel_backend <https://joblib.readthedocs.io/en/latest/parallel.html#joblib.parallel_backend>`_ block.
**kwargs: optional keyword parameters. Any further parameters are passed directly to the distance function. See the docs of the various distance classes in this module.
Returns:
numpy array of shape (nxm): distance matrix
"""
XX = np.reshape(np.arange(len(X)), [-1,1])
- YY = None if Y is None else np.reshape(np.arange(len(Y)), [-1,1])
+ YY = None if Y is None or Y is X else np.reshape(np.arange(len(Y)), [-1,1])
if metric == "bottleneck":
try:
from .. import bottleneck_distance
- return pairwise_distances(XX, YY, metric=_sklearn_wrapper(bottleneck_distance, X, Y, **kwargs))
+ return _pairwise(pairwise_distances, True, XX, YY, metric=_sklearn_wrapper(bottleneck_distance, X, Y, **kwargs), n_jobs=n_jobs)
except ImportError:
print("Gudhi built without CGAL")
raise
elif metric == "pot_wasserstein":
try:
from gudhi.wasserstein import wasserstein_distance as pot_wasserstein_distance
- return pairwise_distances(XX, YY, metric=_sklearn_wrapper(pot_wasserstein_distance, X, Y, **kwargs))
+ return _pairwise(pairwise_distances, True, XX, YY, metric=_sklearn_wrapper(pot_wasserstein_distance, X, Y, **kwargs), n_jobs=n_jobs)
except ImportError:
print("POT (Python Optimal Transport) is not installed. Please install POT or use metric='wasserstein' or metric='hera_wasserstein'")
raise
elif metric == "sliced_wasserstein":
Xproj = _compute_persistence_diagram_projections(X, **kwargs)
Yproj = None if Y is None else _compute_persistence_diagram_projections(Y, **kwargs)
- return pairwise_distances(XX, YY, metric=_sklearn_wrapper(_sliced_wasserstein_distance_on_projections, Xproj, Yproj))
+ return _pairwise(pairwise_distances, True, XX, YY, metric=_sklearn_wrapper(_sliced_wasserstein_distance_on_projections, Xproj, Yproj), n_jobs=n_jobs)
elif type(metric) == str:
- return pairwise_distances(XX, YY, metric=_sklearn_wrapper(PAIRWISE_DISTANCE_FUNCTIONS[metric], X, Y, **kwargs))
+ return _pairwise(pairwise_distances, True, XX, YY, metric=_sklearn_wrapper(PAIRWISE_DISTANCE_FUNCTIONS[metric], X, Y, **kwargs), n_jobs=n_jobs)
else:
- return pairwise_distances(XX, YY, metric=_sklearn_wrapper(metric, X, Y, **kwargs))
+ return _pairwise(pairwise_distances, True, XX, YY, metric=_sklearn_wrapper(metric, X, Y, **kwargs), n_jobs=n_jobs)
class SlicedWassersteinDistance(BaseEstimator, TransformerMixin):
"""
This is a class for computing the sliced Wasserstein distance matrix from a list of persistence diagrams. The Sliced Wasserstein distance is computed by projecting the persistence diagrams onto lines, comparing the projections with the 1-norm, and finally integrating over all possible lines. See http://proceedings.mlr.press/v70/carriere17a.html for more details.
"""
- def __init__(self, num_directions=10):
+ def __init__(self, num_directions=10, n_jobs=None):
"""
Constructor for the SlicedWassersteinDistance class.
Parameters:
num_directions (int): number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the distance computation (default 10).
+ n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_distances` for details.
"""
self.num_directions = num_directions
+ self.n_jobs = n_jobs
def fit(self, X, y=None):
"""
@@ -206,7 +224,7 @@ class SlicedWassersteinDistance(BaseEstimator, TransformerMixin):
Returns:
numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise sliced Wasserstein distances.
"""
- return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="sliced_wasserstein", num_directions=self.num_directions)
+ return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="sliced_wasserstein", num_directions=self.num_directions, n_jobs=self.n_jobs)
def __call__(self, diag1, diag2):
"""
@@ -227,14 +245,16 @@ class BottleneckDistance(BaseEstimator, TransformerMixin):
:Requires: `CGAL <installation.html#cgal>`_ :math:`\geq` 4.11.0
"""
- def __init__(self, epsilon=None):
+ def __init__(self, epsilon=None, n_jobs=None):
"""
Constructor for the BottleneckDistance class.
Parameters:
epsilon (double): absolute (additive) error tolerated on the distance (default is the smallest positive float), see :func:`gudhi.bottleneck_distance`.
+ n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_distances` for details.
"""
self.epsilon = epsilon
+ self.n_jobs = n_jobs
def fit(self, X, y=None):
"""
@@ -257,7 +277,7 @@ class BottleneckDistance(BaseEstimator, TransformerMixin):
Returns:
numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise bottleneck distances.
"""
- Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric="bottleneck", e=self.epsilon)
+ Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric="bottleneck", e=self.epsilon, n_jobs=self.n_jobs)
return Xfit
def __call__(self, diag1, diag2):
@@ -282,15 +302,17 @@ class PersistenceFisherDistance(BaseEstimator, TransformerMixin):
"""
This is a class for computing the persistence Fisher distance matrix from a list of persistence diagrams. The persistence Fisher distance is obtained by computing the original Fisher distance between the probability distributions associated to the persistence diagrams given by convolving them with a Gaussian kernel. See http://papers.nips.cc/paper/8205-persistence-fisher-kernel-a-riemannian-manifold-kernel-for-persistence-diagrams for more details.
"""
- def __init__(self, bandwidth=1., kernel_approx=None):
+ def __init__(self, bandwidth=1., kernel_approx=None, n_jobs=None):
"""
Constructor for the PersistenceFisherDistance class.
Parameters:
bandwidth (double): bandwidth of the Gaussian kernel used to turn persistence diagrams into probability distributions (default 1.).
kernel_approx (class): kernel approximation class used to speed up computation (default None). Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
+ n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_distances` for details.
"""
self.bandwidth, self.kernel_approx = bandwidth, kernel_approx
+ self.n_jobs = n_jobs
def fit(self, X, y=None):
"""
@@ -313,7 +335,7 @@ class PersistenceFisherDistance(BaseEstimator, TransformerMixin):
Returns:
numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence Fisher distances.
"""
- return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="persistence_fisher", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)
+ return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="persistence_fisher", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx, n_jobs=self.n_jobs)
def __call__(self, diag1, diag2):
"""
@@ -332,7 +354,7 @@ class WassersteinDistance(BaseEstimator, TransformerMixin):
"""
This is a class for computing the Wasserstein distance matrix from a list of persistence diagrams.
"""
- def __init__(self, order=2, internal_p=2, mode="pot", delta=0.01):
+ def __init__(self, order=2, internal_p=2, mode="pot", delta=0.01, n_jobs=None):
"""
Constructor for the WassersteinDistance class.
@@ -341,10 +363,12 @@ class WassersteinDistance(BaseEstimator, TransformerMixin):
internal_p (int): ground metric on the (upper-half) plane (i.e. norm l_p in R^2), default value is 2 (euclidean norm), see :func:`gudhi.wasserstein.wasserstein_distance`.
mode (str): method for computing Wasserstein distance. Either "pot" or "hera".
delta (float): relative error 1+delta. Used only if mode == "hera".
+ n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_distances` for details.
"""
self.order, self.internal_p, self.mode = order, internal_p, mode
self.metric = "pot_wasserstein" if mode == "pot" else "hera_wasserstein"
self.delta = delta
+ self.n_jobs = n_jobs
def fit(self, X, y=None):
"""
@@ -368,9 +392,9 @@ class WassersteinDistance(BaseEstimator, TransformerMixin):
numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise Wasserstein distances.
"""
if self.metric == "hera_wasserstein":
- Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p, delta=self.delta)
+ Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p, delta=self.delta, n_jobs=self.n_jobs)
else:
- Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p, matching=False)
+ Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p, matching=False, n_jobs=self.n_jobs)
return Xfit
def __call__(self, diag1, diag2):
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-30 20:37:39 +0000