diff options
Diffstat (limited to 'src/python/gudhi')
| -rw-r--r-- | src/python/gudhi/cubical_complex.pyx | 12 | ||||
| -rw-r--r-- | src/python/gudhi/datasets/generators/_points.cc | 121 | ||||
| -rw-r--r-- | src/python/gudhi/datasets/generators/points.cc | 68 | ||||
| -rw-r--r-- | src/python/gudhi/datasets/generators/points.py | 59 | ||||
| -rw-r--r-- | src/python/gudhi/periodic_cubical_complex.pyx | 12 | ||||
| -rw-r--r-- | src/python/gudhi/point_cloud/knn.py | 10 | ||||
| -rw-r--r-- | src/python/gudhi/representations/vector_methods.py | 80 | ||||
| -rw-r--r-- | src/python/gudhi/simplex_tree.pxd | 2 | ||||
| -rw-r--r-- | src/python/gudhi/simplex_tree.pyx | 31 |
9 files changed, 276 insertions, 119 deletions
diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx index 28fbe3af..8e244bb8 100644 --- a/src/python/gudhi/cubical_complex.pyx +++ b/src/python/gudhi/cubical_complex.pyx @@ -35,7 +35,7 @@ cdef extern from "Cubical_complex_interface.h" namespace "Gudhi": 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) nogil - void compute_persistence(int homology_coeff_field, double min_persistence) nogil + void compute_persistence(int homology_coeff_field, double min_persistence) nogil except+ vector[pair[int, pair[double, double]]] get_persistence() nogil vector[vector[int]] cofaces_of_cubical_persistence_pairs() nogil vector[int] betti_numbers() nogil @@ -147,7 +147,7 @@ cdef class CubicalComplex: :func:`persistence` returns. :param homology_coeff_field: The homology coefficient field. Must be a - prime number + prime number. Default value is 11. Max is 46337. :type homology_coeff_field: int. :param min_persistence: The minimum persistence value to take into account (strictly greater than min_persistence). Default value is @@ -169,7 +169,7 @@ cdef class CubicalComplex: """This function computes and returns the persistence of the complex. :param homology_coeff_field: The homology coefficient field. Must be a - prime number + prime number. Default value is 11. Max is 46337. :type homology_coeff_field: int. :param min_persistence: The minimum persistence value to take into account (strictly greater than min_persistence). Default value is @@ -281,4 +281,8 @@ cdef class CubicalComplex: launched first. """ assert self.pcohptr != NULL, "compute_persistence() must be called before persistence_intervals_in_dimension()" - return np.array(self.pcohptr.intervals_in_dimension(dimension)) + piid = np.array(self.pcohptr.intervals_in_dimension(dimension)) + # Workaround https://github.com/GUDHI/gudhi-devel/issues/507 + if len(piid) == 0: + return np.empty(shape = [0, 2]) + return piid diff --git a/src/python/gudhi/datasets/generators/_points.cc b/src/python/gudhi/datasets/generators/_points.cc new file mode 100644 index 00000000..82fea25b --- /dev/null +++ b/src/python/gudhi/datasets/generators/_points.cc @@ -0,0 +1,121 @@ +/* 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): Hind Montassif + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <pybind11/pybind11.h> +#include <pybind11/numpy.h> + +#include <gudhi/random_point_generators.h> +#include <gudhi/Debug_utils.h> + +#include <CGAL/Epick_d.h> + +namespace py = pybind11; + + +typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kern; + +py::array_t<double> generate_points_on_sphere(size_t n_samples, int ambient_dim, double radius, std::string sample) { + + if (sample != "random") { + throw pybind11::value_error("This sample type is not supported"); + } + + py::array_t<double> points({n_samples, (size_t)ambient_dim}); + + py::buffer_info buf = points.request(); + double *ptr = static_cast<double *>(buf.ptr); + + GUDHI_CHECK(n_samples == buf.shape[0], "Py array first dimension not matching n_samples on sphere"); + GUDHI_CHECK(ambient_dim == buf.shape[1], "Py array second dimension not matching the ambient space dimension"); + + + std::vector<typename Kern::Point_d> points_generated; + + { + py::gil_scoped_release release; + points_generated = Gudhi::generate_points_on_sphere_d<Kern>(n_samples, ambient_dim, radius); + } + + for (size_t i = 0; i < n_samples; i++) + for (int j = 0; j < ambient_dim; j++) + ptr[i*ambient_dim+j] = points_generated[i][j]; + + return points; +} + +py::array_t<double> generate_points_on_torus(size_t n_samples, int dim, std::string sample) { + + if ( (sample != "random") && (sample != "grid")) { + throw pybind11::value_error("This sample type is not supported"); + } + + std::vector<typename Kern::Point_d> points_generated; + + { + py::gil_scoped_release release; + points_generated = Gudhi::generate_points_on_torus_d<Kern>(n_samples, dim, sample); + } + + size_t npoints = points_generated.size(); + + GUDHI_CHECK(2*dim == points_generated[0].size(), "Py array second dimension not matching the double torus dimension"); + + py::array_t<double> points({npoints, (size_t)2*dim}); + + py::buffer_info buf = points.request(); + double *ptr = static_cast<double *>(buf.ptr); + + for (size_t i = 0; i < npoints; i++) + for (int j = 0; j < 2*dim; j++) + ptr[i*(2*dim)+j] = points_generated[i][j]; + + return points; +} + +PYBIND11_MODULE(_points, m) { + m.attr("__license__") = "LGPL v3"; + + m.def("sphere", &generate_points_on_sphere, + py::arg("n_samples"), py::arg("ambient_dim"), + py::arg("radius") = 1., py::arg("sample") = "random", + R"pbdoc( + Generate random i.i.d. points uniformly on a (d-1)-sphere in R^d + + :param n_samples: The number of points to be generated. + :type n_samples: integer + :param ambient_dim: The ambient dimension d. + :type ambient_dim: integer + :param radius: The radius. Default value is `1.`. + :type radius: float + :param sample: The sample type. Default and only available value is `"random"`. + :type sample: string + :returns: the generated points on a sphere. + )pbdoc"); + + m.def("ctorus", &generate_points_on_torus, + py::arg("n_samples"), py::arg("dim"), py::arg("sample") = "random", + R"pbdoc( + Generate random i.i.d. points on a d-torus in R^2d or as a grid + + :param n_samples: The number of points to be generated. + :type n_samples: integer + :param dim: The dimension of the torus on which points would be generated in R^2*dim. + :type dim: integer + :param sample: The sample type. Available values are: `"random"` and `"grid"`. Default value is `"random"`. + :type sample: string + :returns: the generated points on a torus. + + The shape of returned numpy array is: + + If sample is 'random': (n_samples, 2*dim). + + If sample is 'grid': (⌊n_samples**(1./dim)⌋**dim, 2*dim), where shape[0] is rounded down to the closest perfect 'dim'th power. + )pbdoc"); +} diff --git a/src/python/gudhi/datasets/generators/points.cc b/src/python/gudhi/datasets/generators/points.cc deleted file mode 100644 index d658946b..00000000 --- a/src/python/gudhi/datasets/generators/points.cc +++ /dev/null @@ -1,68 +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): Hind Montassif - * - * Copyright (C) 2021 Inria - * - * Modification(s): - * - YYYY/MM Author: Description of the modification - */ - -#include <pybind11/pybind11.h> -#include <pybind11/numpy.h> - -#include <gudhi/random_point_generators.h> -#include <gudhi/Debug_utils.h> - -#include <CGAL/Epick_d.h> - -namespace py = pybind11; - - -typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kern; - -py::array_t<double> generate_points_on_sphere(size_t n_samples, int ambient_dim, double radius, std::string sample) { - - if (sample != "random") { - throw pybind11::value_error("This sample type is not supported"); - } - - py::array_t<double> points({n_samples, (size_t)ambient_dim}); - - py::buffer_info buf = points.request(); - double *ptr = static_cast<double *>(buf.ptr); - - GUDHI_CHECK(n_samples == buf.shape[0], "Py array first dimension not matching n_samples on sphere"); - GUDHI_CHECK(ambient_dim == buf.shape[1], "Py array second dimension not matching the ambient space dimension"); - - - py::gil_scoped_release release; - auto points_generated = Gudhi::generate_points_on_sphere_d<Kern>(n_samples, ambient_dim, radius); - - for (size_t i = 0; i < n_samples; i++) - for (int j = 0; j < ambient_dim; j++) - ptr[i*ambient_dim+j] = points_generated[i][j]; - - return points; -} - -PYBIND11_MODULE(points, m) { - m.attr("__license__") = "LGPL v3"; - m.def("sphere", &generate_points_on_sphere, - py::arg("n_samples"), py::arg("ambient_dim"), - py::arg("radius") = 1., py::arg("sample") = "random", - R"pbdoc( - Generate random i.i.d. points uniformly on a (d-1)-sphere in R^d - - :param n_samples: The number of points to be generated. - :type n_samples: integer - :param ambient_dim: The ambient dimension d. - :type ambient_dim: integer - :param radius: The radius. Default value is `1.`. - :type radius: float - :param sample: The sample type. Default and only available value is `"random"`. - :type sample: string - :rtype: numpy array of float - :returns: the generated points on a sphere. - )pbdoc"); -} diff --git a/src/python/gudhi/datasets/generators/points.py b/src/python/gudhi/datasets/generators/points.py new file mode 100644 index 00000000..9bb2799d --- /dev/null +++ b/src/python/gudhi/datasets/generators/points.py @@ -0,0 +1,59 @@ +# 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): Hind Montassif +# +# Copyright (C) 2021 Inria +# +# Modification(s): +# - YYYY/MM Author: Description of the modification + +import numpy as np + +from ._points import ctorus +from ._points import sphere + +def _generate_random_points_on_torus(n_samples, dim): + + # Generate random angles of size n_samples*dim + alpha = 2*np.pi*np.random.rand(n_samples*dim) + + # Based on angles, construct points of size n_samples*dim on a circle and reshape the result in a n_samples*2*dim array + array_points = np.column_stack([np.cos(alpha), np.sin(alpha)]).reshape(-1, 2*dim) + + return array_points + +def _generate_grid_points_on_torus(n_samples, dim): + + # Generate points on a dim-torus as a grid + n_samples_grid = int((n_samples+.5)**(1./dim)) # add .5 to avoid rounding down with numerical approximations + alpha = np.linspace(0, 2*np.pi, n_samples_grid, endpoint=False) + + array_points = np.column_stack([np.cos(alpha), np.sin(alpha)]) + array_points_idx = np.empty([n_samples_grid]*dim + [dim], dtype=int) + for i, x in enumerate(np.ix_(*([np.arange(n_samples_grid)]*dim))): + array_points_idx[...,i] = x + return array_points[array_points_idx].reshape(-1, 2*dim) + +def torus(n_samples, dim, sample='random'): + """ + Generate points on a flat dim-torus in R^2dim either randomly or on a grid + + :param n_samples: The number of points to be generated. + :param dim: The dimension of the torus on which points would be generated in R^2*dim. + :param sample: The sample type of the generated points. Can be 'random' or 'grid'. + :returns: numpy array containing the generated points on a torus. + + The shape of returned numpy array is: + + If sample is 'random': (n_samples, 2*dim). + + If sample is 'grid': (⌊n_samples**(1./dim)⌋**dim, 2*dim), where shape[0] is rounded down to the closest perfect 'dim'th power. + """ + if sample == 'random': + # Generate points randomly + return _generate_random_points_on_torus(n_samples, dim) + elif sample == 'grid': + # Generate points on a grid + return _generate_grid_points_on_torus(n_samples, dim) + else: + raise ValueError("Sample type '{}' is not supported".format(sample)) diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx index d353d2af..6c21e902 100644 --- a/src/python/gudhi/periodic_cubical_complex.pyx +++ b/src/python/gudhi/periodic_cubical_complex.pyx @@ -32,7 +32,7 @@ cdef extern from "Cubical_complex_interface.h" namespace "Gudhi": 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) nogil - void compute_persistence(int homology_coeff_field, double min_persistence) nogil + void compute_persistence(int homology_coeff_field, double min_persistence) nogil except + vector[pair[int, pair[double, double]]] get_persistence() nogil vector[vector[int]] cofaces_of_cubical_persistence_pairs() nogil vector[int] betti_numbers() nogil @@ -148,7 +148,7 @@ cdef class PeriodicCubicalComplex: :func:`persistence` returns. :param homology_coeff_field: The homology coefficient field. Must be a - prime number + prime number. Default value is 11. Max is 46337. :type homology_coeff_field: int. :param min_persistence: The minimum persistence value to take into account (strictly greater than min_persistence). Default value is @@ -170,7 +170,7 @@ cdef class PeriodicCubicalComplex: """This function computes and returns the persistence of the complex. :param homology_coeff_field: The homology coefficient field. Must be a - prime number + prime number. Default value is 11. Max is 46337. :type homology_coeff_field: int. :param min_persistence: The minimum persistence value to take into account (strictly greater than min_persistence). Default value is @@ -280,4 +280,8 @@ cdef class PeriodicCubicalComplex: launched first. """ assert self.pcohptr != NULL, "compute_persistence() must be called before persistence_intervals_in_dimension()" - return np.array(self.pcohptr.intervals_in_dimension(dimension)) + piid = np.array(self.pcohptr.intervals_in_dimension(dimension)) + # Workaround https://github.com/GUDHI/gudhi-devel/issues/507 + if len(piid) == 0: + return np.empty(shape = [0, 2]) + return piid diff --git a/src/python/gudhi/point_cloud/knn.py b/src/python/gudhi/point_cloud/knn.py index 829bf1bf..de5844f9 100644 --- a/src/python/gudhi/point_cloud/knn.py +++ b/src/python/gudhi/point_cloud/knn.py @@ -8,6 +8,7 @@ # - YYYY/MM Author: Description of the modification import numpy +import warnings # TODO: https://github.com/facebookresearch/faiss @@ -257,6 +258,9 @@ class KNearestNeighbors: if ef is not None: self.graph.set_ef(ef) neighbors, distances = self.graph.knn_query(X, k, num_threads=self.params["num_threads"]) + with warnings.catch_warnings(): + if not(numpy.all(numpy.isfinite(distances))): + warnings.warn("Overflow/infinite value encountered while computing 'distances'", RuntimeWarning) # The k nearest neighbors are always sorted. I couldn't find it in the doc, but the code calls searchKnn, # which returns a priority_queue, and then fills the return array backwards with top/pop on the queue. if self.return_index: @@ -290,6 +294,9 @@ class KNearestNeighbors: if self.return_index: if self.return_distance: distances, neighbors = mat.Kmin_argKmin(k, dim=1) + with warnings.catch_warnings(): + if not(torch.isfinite(distances).all()): + warnings.warn("Overflow/infinite value encountered while computing 'distances'", RuntimeWarning) if p != numpy.inf: distances = distances ** (1.0 / p) return neighbors, distances @@ -298,6 +305,9 @@ class KNearestNeighbors: return neighbors if self.return_distance: distances = mat.Kmin(k, dim=1) + with warnings.catch_warnings(): + if not(torch.isfinite(distances).all()): + warnings.warn("Overflow/infinite value encountered while computing 'distances'", RuntimeWarning) if p != numpy.inf: distances = distances ** (1.0 / p) return distances diff --git a/src/python/gudhi/representations/vector_methods.py b/src/python/gudhi/representations/vector_methods.py index 84bc99a2..e883b5dd 100644 --- a/src/python/gudhi/representations/vector_methods.py +++ b/src/python/gudhi/representations/vector_methods.py @@ -6,6 +6,7 @@ # # Modification(s): # - 2020/06 Martin: ATOL integration +# - 2021/11 Vincent Rouvreau: factorize _automatic_sample_range import numpy as np from sklearn.base import BaseEstimator, TransformerMixin @@ -45,10 +46,14 @@ class PersistenceImage(BaseEstimator, TransformerMixin): y (n x 1 array): persistence diagram labels (unused). """ if np.isnan(np.array(self.im_range)).any(): - new_X = BirthPersistenceTransform().fit_transform(X) - pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(new_X,y) - [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]] - self.im_range = np.where(np.isnan(np.array(self.im_range)), np.array([mx, Mx, my, My]), np.array(self.im_range)) + try: + new_X = BirthPersistenceTransform().fit_transform(X) + pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(new_X,y) + [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]] + self.im_range = np.where(np.isnan(np.array(self.im_range)), np.array([mx, Mx, my, My]), np.array(self.im_range)) + except ValueError: + # Empty persistence diagram case - https://github.com/GUDHI/gudhi-devel/issues/507 + pass return self def transform(self, X): @@ -94,6 +99,28 @@ class PersistenceImage(BaseEstimator, TransformerMixin): """ return self.fit_transform([diag])[0,:] +def _automatic_sample_range(sample_range, X, y): + """ + Compute and returns sample range from the persistence diagrams if one of the sample_range values is numpy.nan. + + Parameters: + sample_range (a numpy array of 2 float): minimum and maximum of all piecewise-linear function domains, of + the form [x_min, x_max]. + X (list of n x 2 numpy arrays): input persistence diagrams. + y (n x 1 array): persistence diagram labels (unused). + """ + nan_in_range = np.isnan(sample_range) + if nan_in_range.any(): + try: + pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(X,y) + [mx,my] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]] + [Mx,My] = [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]] + return np.where(nan_in_range, np.array([mx, My]), sample_range) + except ValueError: + # Empty persistence diagram case - https://github.com/GUDHI/gudhi-devel/issues/507 + pass + return sample_range + class Landscape(BaseEstimator, TransformerMixin): """ This is a class for computing persistence landscapes from a list of persistence diagrams. A persistence landscape is a collection of 1D piecewise-linear functions computed from the rank function associated to the persistence diagram. These piecewise-linear functions are then sampled evenly on a given range and the corresponding vectors of samples are concatenated and returned. See http://jmlr.org/papers/v16/bubenik15a.html for more details. @@ -119,10 +146,7 @@ class Landscape(BaseEstimator, TransformerMixin): X (list of n x 2 numpy arrays): input persistence diagrams. y (n x 1 array): persistence diagram labels (unused). """ - if self.nan_in_range.any(): - pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(X,y) - [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]] - self.sample_range = np.where(self.nan_in_range, np.array([mx, My]), np.array(self.sample_range)) + self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y) return self def transform(self, X): @@ -218,10 +242,7 @@ class Silhouette(BaseEstimator, TransformerMixin): X (list of n x 2 numpy arrays): input persistence diagrams. y (n x 1 array): persistence diagram labels (unused). """ - if np.isnan(np.array(self.sample_range)).any(): - pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(X,y) - [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]] - self.sample_range = np.where(np.isnan(np.array(self.sample_range)), np.array([mx, My]), np.array(self.sample_range)) + self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y) return self def transform(self, X): @@ -307,10 +328,7 @@ class BettiCurve(BaseEstimator, TransformerMixin): X (list of n x 2 numpy arrays): input persistence diagrams. y (n x 1 array): persistence diagram labels (unused). """ - if np.isnan(np.array(self.sample_range)).any(): - pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(X,y) - [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]] - self.sample_range = np.where(np.isnan(np.array(self.sample_range)), np.array([mx, My]), np.array(self.sample_range)) + self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y) return self def transform(self, X): @@ -374,10 +392,7 @@ class Entropy(BaseEstimator, TransformerMixin): X (list of n x 2 numpy arrays): input persistence diagrams. y (n x 1 array): persistence diagram labels (unused). """ - if np.isnan(np.array(self.sample_range)).any(): - pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(X,y) - [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]] - self.sample_range = np.where(np.isnan(np.array(self.sample_range)), np.array([mx, My]), np.array(self.sample_range)) + self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y) return self def transform(self, X): @@ -396,9 +411,13 @@ class Entropy(BaseEstimator, TransformerMixin): new_X = BirthPersistenceTransform().fit_transform(X) for i in range(num_diag): - orig_diagram, diagram, num_pts_in_diag = X[i], new_X[i], X[i].shape[0] - new_diagram = DiagramScaler(use=True, scalers=[([1], MaxAbsScaler())]).fit_transform([diagram])[0] + try: + new_diagram = DiagramScaler(use=True, scalers=[([1], MaxAbsScaler())]).fit_transform([diagram])[0] + except ValueError: + # Empty persistence diagram case - https://github.com/GUDHI/gudhi-devel/issues/507 + assert len(diagram) == 0 + new_diagram = np.empty(shape = [0, 2]) if self.mode == "scalar": ent = - np.sum( np.multiply(new_diagram[:,1], np.log(new_diagram[:,1])) ) @@ -412,12 +431,11 @@ class Entropy(BaseEstimator, TransformerMixin): max_idx = np.clip(np.ceil((py - self.sample_range[0]) / step_x).astype(int), 0, self.resolution) for k in range(min_idx, max_idx): ent[k] += (-1) * new_diagram[j,1] * np.log(new_diagram[j,1]) - if self.normalized: - ent = ent / np.linalg.norm(ent, ord=1) - Xfit.append(np.reshape(ent,[1,-1])) - - Xfit = np.concatenate(Xfit, 0) + if self.normalized: + ent = ent / np.linalg.norm(ent, ord=1) + Xfit.append(np.reshape(ent,[1,-1])) + Xfit = np.concatenate(Xfit, axis=0) return Xfit def __call__(self, diag): @@ -478,7 +496,13 @@ class TopologicalVector(BaseEstimator, TransformerMixin): diagram, num_pts_in_diag = X[i], X[i].shape[0] pers = 0.5 * (diagram[:,1]-diagram[:,0]) min_pers = np.minimum(pers,np.transpose(pers)) - distances = DistanceMetric.get_metric("chebyshev").pairwise(diagram) + # Works fine with sklearn 1.0, but an ValueError exception is thrown on past versions + try: + distances = DistanceMetric.get_metric("chebyshev").pairwise(diagram) + except ValueError: + # Empty persistence diagram case - https://github.com/GUDHI/gudhi-devel/issues/507 + assert len(diagram) == 0 + distances = np.empty(shape = [0, 0]) vect = np.flip(np.sort(np.triu(np.minimum(distances, min_pers)), axis=None), 0) dim = min(len(vect), thresh) Xfit[i, :dim] = vect[:dim] diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 3b8ea4f9..006a24ed 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -78,7 +78,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>": Simplex_tree_persistence_interface(Simplex_tree_interface_full_featured * st, bool persistence_dim_max) nogil - void compute_persistence(int homology_coeff_field, double min_persistence) nogil + void compute_persistence(int homology_coeff_field, double min_persistence) nogil except + vector[pair[int, pair[double, double]]] get_persistence() nogil vector[int] betti_numbers() nogil vector[int] persistent_betti_numbers(double from_value, double to_value) nogil diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index be08a3a1..c3720936 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -9,8 +9,7 @@ from cython.operator import dereference, preincrement from libc.stdint cimport intptr_t -import numpy -from numpy import array as np_array +import numpy as np cimport gudhi.simplex_tree __author__ = "Vincent Rouvreau" @@ -412,7 +411,7 @@ cdef class SimplexTree: """This function retrieves good values for extended persistence, and separate the diagrams into the Ordinary, Relative, Extended+ and Extended- subdiagrams. - :param homology_coeff_field: The homology coefficient field. Must be a prime number. Default value is 11. + :param homology_coeff_field: The homology coefficient field. Must be a prime number. Default value is 11. Max is 46337. :type homology_coeff_field: int :param min_persistence: The minimum persistence value (i.e., the absolute value of the difference between the persistence diagram point coordinates) to take into account (strictly greater than min_persistence). @@ -449,7 +448,7 @@ cdef class SimplexTree: """This function computes and returns the persistence of the simplicial complex. :param homology_coeff_field: The homology coefficient field. Must be a - prime number. Default value is 11. + prime number. Default value is 11. Max is 46337. :type homology_coeff_field: int :param min_persistence: The minimum persistence value to take into account (strictly greater than min_persistence). Default value is @@ -472,7 +471,7 @@ cdef class SimplexTree: when you do not want the list :func:`persistence` returns. :param homology_coeff_field: The homology coefficient field. Must be a - prime number. Default value is 11. + prime number. Default value is 11. Max is 46337. :type homology_coeff_field: int :param min_persistence: The minimum persistence value to take into account (strictly greater than min_persistence). Default value is @@ -542,7 +541,11 @@ cdef class SimplexTree: function to be launched first. """ assert self.pcohptr != NULL, "compute_persistence() must be called before persistence_intervals_in_dimension()" - return np_array(self.pcohptr.intervals_in_dimension(dimension)) + piid = np.array(self.pcohptr.intervals_in_dimension(dimension)) + # Workaround https://github.com/GUDHI/gudhi-devel/issues/507 + if len(piid) == 0: + return np.empty(shape = [0, 2]) + return piid def persistence_pairs(self): """This function returns a list of persistence birth and death simplices pairs. @@ -583,8 +586,8 @@ cdef class SimplexTree: """ assert self.pcohptr != NULL, "lower_star_persistence_generators() requires that persistence() be called first." gen = self.pcohptr.lower_star_generators() - normal = [np_array(d).reshape(-1,2) for d in gen.first] - infinite = [np_array(d) for d in gen.second] + normal = [np.array(d).reshape(-1,2) for d in gen.first] + infinite = [np.array(d) for d in gen.second] return (normal, infinite) def flag_persistence_generators(self): @@ -602,19 +605,19 @@ cdef class SimplexTree: assert self.pcohptr != NULL, "flag_persistence_generators() requires that persistence() be called first." gen = self.pcohptr.flag_generators() if len(gen.first) == 0: - normal0 = numpy.empty((0,3)) + normal0 = np.empty((0,3)) normals = [] else: l = iter(gen.first) - normal0 = np_array(next(l)).reshape(-1,3) - normals = [np_array(d).reshape(-1,4) for d in l] + normal0 = np.array(next(l)).reshape(-1,3) + normals = [np.array(d).reshape(-1,4) for d in l] if len(gen.second) == 0: - infinite0 = numpy.empty(0) + infinite0 = np.empty(0) infinites = [] else: l = iter(gen.second) - infinite0 = np_array(next(l)) - infinites = [np_array(d).reshape(-1,2) for d in l] + infinite0 = np.array(next(l)) + infinites = [np.array(d).reshape(-1,2) for d in l] return (normal0, normals, infinite0, infinites) def collapse_edges(self, nb_iterations = 1): |