diff options
| -rw-r--r-- | .github/for_maintainers/new_gudhi_version_creation.md | 10 | ||||
| -rw-r--r-- | .github/next_release.md | 50 | ||||
| -rw-r--r-- | CMakeGUDHIVersion.txt | 4 | ||||
| -rw-r--r-- | Dockerfile_gudhi_installation | 6 | ||||
| -rw-r--r-- | src/cmake/modules/GUDHI_user_version_target.cmake | 13 | ||||
| -rw-r--r-- | src/python/doc/point_cloud_sum.inc | 4 | ||||
| -rw-r--r-- | src/python/gudhi/cubical_complex.pyx | 50 | ||||
| -rw-r--r-- | src/python/gudhi/periodic_cubical_complex.pyx | 86 | ||||
| -rw-r--r-- | src/python/gudhi/point_cloud/dtm.py | 109 | ||||
| -rw-r--r-- | src/python/gudhi/point_cloud/knn.py | 4 | ||||
| -rwxr-xr-x | src/python/test/test_cubical_complex.py | 17 | ||||
| -rwxr-xr-x | src/python/test/test_dtm.py | 23 |
12 files changed, 262 insertions, 114 deletions
diff --git a/.github/for_maintainers/new_gudhi_version_creation.md b/.github/for_maintainers/new_gudhi_version_creation.md index 8674222b..86c393a0 100644 --- a/.github/for_maintainers/new_gudhi_version_creation.md +++ b/.github/for_maintainers/new_gudhi_version_creation.md @@ -17,8 +17,7 @@ rm -rf data/points/COIL_database/lucky_cat.off_dist data/points/COIL_database/lu Checkin the modifications, build and test the version: ```bash git submodule update --init -mkdir build -cd build +rm -rf build; mkdir build; cd build cmake -DCGAL_DIR=/your/path/to/CGAL -DWITH_GUDHI_EXAMPLE=ON -DWITH_GUDHI_BENCHMARK=ON -DUSER_VERSION_DIR=gudhi.@GUDHI_VERSION@ -DPython_ADDITIONAL_VERSIONS=3 .. make user_version date +"%d-%m-%Y-%T" > gudhi.@GUDHI_VERSION@/timestamp.txt @@ -27,7 +26,7 @@ md5sum gudhi.@GUDHI_VERSION@.tar.gz > md5sum.txt sha256sum gudhi.@GUDHI_VERSION@.tar.gz > sha256sum.txt sha512sum gudhi.@GUDHI_VERSION@.tar.gz > sha512sum.txt -make -j all test +make -j 4 all && ctest -j 4 --output-on-failure ``` ***[Check there are no error]*** @@ -43,7 +42,8 @@ make doxygen 2>&1 | tee dox.log && grep warning dox.log ```bash cp -R gudhi.@GUDHI_VERSION@/doc/html gudhi.doc.@GUDHI_VERSION@/cpp cd gudhi.@GUDHI_VERSION@ -rm -rf build; mkdir build; cd build; cmake -DCGAL_DIR=/your/path/to/CGAL -DWITH_GUDHI_EXAMPLE=ON -DPython_ADDITIONAL_VERSIONS=3 .. +rm -rf build; mkdir build; cd build +cmake -DCGAL_DIR=/your/path/to/CGAL -DWITH_GUDHI_EXAMPLE=ON -DPython_ADDITIONAL_VERSIONS=3 .. export LC_ALL=en_US.UTF-8 # cf. bug make sphinx ``` @@ -56,7 +56,7 @@ cd ../.. tar -czvf gudhi.doc.@GUDHI_VERSION@.tar.gz gudhi.doc.@GUDHI_VERSION@ cd gudhi.@GUDHI_VERSION@/build -make all test +make -j 4 all && ctest -j 4 --output-on-failure ``` ***[Check there are no error]*** diff --git a/.github/next_release.md b/.github/next_release.md index d3c9ce68..a2805a55 100644 --- a/.github/next_release.md +++ b/.github/next_release.md @@ -1,50 +1,19 @@ -We are pleased to announce the release 3.2.0 of the GUDHI library. +We are pleased to announce the release 3.X.X of the GUDHI library. -As a major new feature, the GUDHI library now offers a Python interface to [Hera](https://bitbucket.org/grey_narn/hera/src/master/) to compute the Wasserstein distance. -[PyBind11](https://github.com/pybind/pybind11) is now required to build the Python module. +As a major new feature, the GUDHI library now offers ... -We are now using GitHub to develop the GUDHI library, do not hesitate to [fork the GUDHI project on GitHub](https://github.com/GUDHI/gudhi-devel). From a user point of view, we recommend to download GUDHI user version (gudhi.3.2.0.tar.gz). +We are now using GitHub to develop the GUDHI library, do not hesitate to [fork the GUDHI project on GitHub](https://github.com/GUDHI/gudhi-devel). From a user point of view, we recommend to download GUDHI user version (gudhi.3.X.X.tar.gz). -Below is a list of changes made since GUDHI 3.1.1: +Below is a list of changes made since GUDHI 3.X-1.X-1: -- Point cloud utilities - - A new module [Time Delay Embedding](https://gudhi.inria.fr/python/latest/point_cloud.html#time-delay-embedding) - to embed time-series data in the R^d according to [Takens' Embedding Theorem](https://en.wikipedia.org/wiki/Takens%27s_theorem) - and obtain the coordinates of each point. - - A new module [K Nearest Neighbors](https://gudhi.inria.fr/python/latest/point_cloud.html#k-nearest-neighbors) - that wraps several implementations for computing the k nearest neighbors in a point set. - - A new module [Distance To Measure](https://gudhi.inria.fr/python/latest/point_cloud.html#distance-to-measure) - to compute the distance to the empirical measure defined by a point set +- [Module](link) + - ... -- [Persistence representations](https://gudhi.inria.fr/python/latest/representations.html) - - Interface to Wasserstein distances. - -- Rips complex - - A new module [Weighted Rips Complex](https://gudhi.inria.fr/python/latest/rips_complex_user.html#weighted-rips-complex) - to construct a simplicial complex from a distance matrix and weights on vertices. - -- [Wassertein distance](https://gudhi.inria.fr/python/latest/wasserstein_distance_user.html) - - An [another implementation](https://gudhi.inria.fr/python/latest/wasserstein_distance_user.html#hera) - comes from Hera (BSD-3-Clause) which is based on [Geometry Helps to Compare Persistence Diagrams](http://doi.acm.org/10.1145/3064175) - by Michael Kerber, Dmitriy Morozov, and Arnur Nigmetov. - - `gudhi.wasserstein.wasserstein_distance` has now an option to return the optimal matching that achieves the distance between the two diagrams. - - A new module [Barycenters](https://gudhi.inria.fr/python/latest/wasserstein_distance_user.html#barycenters) - to estimate the Frechet mean (aka Wasserstein barycenter) between persistence diagrams. - -- [Simplex tree](https://gudhi.inria.fr/python/latest/simplex_tree_ref.html) - - Extend filtration method to compute extended persistence - - Flag and lower star persistence pairs generators - - A new interface to filtration, simplices and skeleton getters to return an iterator - -- [Alpha complex](https://gudhi.inria.fr/doc/latest/group__alpha__complex.html) - - Improve computations (cache circumcenters computation and point comparison improvement) - -- [Persistence graphical tools](https://gudhi.inria.fr/python/latest/persistence_graphical_tools_user.html) - - New rendering option proposed (use LaTeX style, add grey block, improved positioning of labels, etc.). - - Can now handle (N x 2) numpy arrays as input +- [Module](link) + - ... - Miscellaneous - - The [list of bugs that were solved since GUDHI-3.2.0](https://github.com/GUDHI/gudhi-devel/issues?q=label%3A3.2.0+is%3Aclosed) is available on GitHub. + - The [list of bugs that were solved since GUDHI-3.X-1.X-1](https://github.com/GUDHI/gudhi-devel/issues?q=label%3A3.1.1+is%3Aclosed) is available on GitHub. All modules are distributed under the terms of the MIT license. However, there are still GPL dependencies for many modules. We invite you to check our [license dedicated web page](https://gudhi.inria.fr/licensing/) for further details. @@ -56,3 +25,4 @@ We provide [bibtex entries](https://gudhi.inria.fr/doc/latest/_citation.html) fo Feel free to [contact us](https://gudhi.inria.fr/contact/) in case you have any questions or remarks. For further information about downloading and installing the library ([C++](https://gudhi.inria.fr/doc/latest/installation.html) or [Python](https://gudhi.inria.fr/python/latest/installation.html)), please visit the [GUDHI web site](https://gudhi.inria.fr/). + diff --git a/CMakeGUDHIVersion.txt b/CMakeGUDHIVersion.txt index ac89fa4d..b9895bb1 100644 --- a/CMakeGUDHIVersion.txt +++ b/CMakeGUDHIVersion.txt @@ -1,6 +1,6 @@ set (GUDHI_MAJOR_VERSION 3) -set (GUDHI_MINOR_VERSION 2) -set (GUDHI_PATCH_VERSION 0) +set (GUDHI_MINOR_VERSION 3) +set (GUDHI_PATCH_VERSION 0.rc1) set(GUDHI_VERSION ${GUDHI_MAJOR_VERSION}.${GUDHI_MINOR_VERSION}.${GUDHI_PATCH_VERSION}) message(STATUS "GUDHI version : ${GUDHI_VERSION}") diff --git a/Dockerfile_gudhi_installation b/Dockerfile_gudhi_installation index f9e8813b..461a8a19 100644 --- a/Dockerfile_gudhi_installation +++ b/Dockerfile_gudhi_installation @@ -58,9 +58,9 @@ RUN pip3 install \ # apt clean up RUN apt autoremove && rm -rf /var/lib/apt/lists/* -RUN curl -LO "https://github.com/GUDHI/gudhi-devel/releases/download/tags%2Fgudhi-release-3.1.1/gudhi.3.1.1.tar.gz" \ -&& tar xf gudhi.3.1.1.tar.gz \ -&& cd gudhi.3.1.1 \ +RUN curl -LO "https://github.com/GUDHI/gudhi-devel/releases/download/tags%2Fgudhi-release-3.2.0/gudhi.3.2.0.tar.gz" \ +&& tar xf gudhi.3.2.0.tar.gz \ +&& cd gudhi.3.2.0 \ && mkdir build && cd build && cmake -DCMAKE_BUILD_TYPE=Release -DWITH_GUDHI_PYTHON=OFF -DPython_ADDITIONAL_VERSIONS=3 .. \ && make all test install \ && cmake -DWITH_GUDHI_PYTHON=ON . \ diff --git a/src/cmake/modules/GUDHI_user_version_target.cmake b/src/cmake/modules/GUDHI_user_version_target.cmake index 9cf648e3..e99bb42d 100644 --- a/src/cmake/modules/GUDHI_user_version_target.cmake +++ b/src/cmake/modules/GUDHI_user_version_target.cmake @@ -49,8 +49,17 @@ add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_SOURCE_DIR}/CMakeGUDHIVersion.txt ${GUDHI_USER_VERSION_DIR}/CMakeGUDHIVersion.txt) -add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E - copy_directory ${CMAKE_SOURCE_DIR}/${GUDHI_PYTHON_PATH} ${GUDHI_USER_VERSION_DIR}/python) +# As cython generates .cpp files in source, we have to copy all except cpp files from python directory +file(GLOB_RECURSE PYTHON_FILES ${CMAKE_SOURCE_DIR}/${GUDHI_PYTHON_PATH}/*) +foreach(PYTHON_FILE ${PYTHON_FILES}) + get_filename_component(PYTHON_FILE_EXT ${PYTHON_FILE} EXT) + if (NOT "${PYTHON_FILE_EXT}" STREQUAL ".cpp") + string(REPLACE "${CMAKE_SOURCE_DIR}/${GUDHI_PYTHON_PATH}/" "" RELATIVE_PYTHON_FILE ${PYTHON_FILE}) + add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E + copy ${PYTHON_FILE} ${GUDHI_USER_VERSION_DIR}/python/${RELATIVE_PYTHON_FILE}) + endif() +endforeach() + add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_SOURCE_DIR}/data ${GUDHI_USER_VERSION_DIR}/data) add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E diff --git a/src/python/doc/point_cloud_sum.inc b/src/python/doc/point_cloud_sum.inc index 4315cea6..f955c3ab 100644 --- a/src/python/doc/point_cloud_sum.inc +++ b/src/python/doc/point_cloud_sum.inc @@ -3,8 +3,8 @@ +-----------------------------------+---------------------------------------------------------------+-------------------------------------------------------------------+ | | :math:`(x_1, x_2, \ldots, x_d)` | Utilities to process point clouds: read from file, subsample, | :Authors: Vincent Rouvreau, Marc Glisse, Masatoshi Takenouchi | - | | :math:`(y_1, y_2, \ldots, y_d)` | find neighbors, embed time series in higher dimension, etc. | | - | | | :Since: GUDHI 2.0.0 | + | | :math:`(y_1, y_2, \ldots, y_d)` | find neighbors, embed time series in higher dimension, | | + | | estimate a density, etc. | :Since: GUDHI 2.0.0 | | | | | | | | :License: MIT (`GPL v3 </licensing/>`_, BSD-3-Clause, Apache-2.0) | +-----------------------------------+---------------------------------------------------------------+-------------------------------------------------------------------+ 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/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/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py index 5c59db8f..d0e4e9e8 100755 --- a/src/python/test/test_cubical_complex.py +++ b/src/python/test/test_cubical_complex.py @@ -157,3 +157,20 @@ def test_cubical_generators(): assert np.array_equal(g[0][0], np.empty(shape=[0,2])) assert np.array_equal(g[0][1], np.array([[7, 4]])) assert np.array_equal(g[1][0], np.array([8])) + +def test_cubical_cofaces_of_persistence_pairs_when_pd_has_no_paired_birth_and_death(): + cubCpx = CubicalComplex(dimensions=[1,2], top_dimensional_cells=[0.0, 1.0]) + Diag = cubCpx.persistence(homology_coeff_field=2, min_persistence=0) + pairs = cubCpx.cofaces_of_persistence_pairs() + assert pairs[0] == [] + assert np.array_equal(pairs[1][0], np.array([0])) + +def test_periodic_cofaces_of_persistence_pairs_when_pd_has_no_paired_birth_and_death(): + perCubCpx = PeriodicCubicalComplex(dimensions=[1,2], top_dimensional_cells=[0.0, 1.0], + periodic_dimensions=[True, True]) + Diag = perCubCpx.persistence(homology_coeff_field=2, min_persistence=0) + pairs = perCubCpx.cofaces_of_persistence_pairs() + assert pairs[0] == [] + assert np.array_equal(pairs[1][0], np.array([0])) + assert np.array_equal(pairs[1][1], np.array([0, 1])) + assert np.array_equal(pairs[1][2], np.array([1])) diff --git a/src/python/test/test_dtm.py b/src/python/test/test_dtm.py index bff4c267..0a52279e 100755 --- a/src/python/test/test_dtm.py +++ b/src/python/test/test_dtm.py @@ -8,10 +8,11 @@ - YYYY/MM Author: Description of the modification """ -from gudhi.point_cloud.dtm import DistanceToMeasure +from gudhi.point_cloud.dtm import DistanceToMeasure, DTMDensity import numpy import pytest import torch +import math def test_dtm_compare_euclidean(): @@ -66,3 +67,23 @@ def test_dtm_precomputed(): dtm = DistanceToMeasure(2, q=2, metric="neighbors") r = dtm.fit_transform(dist) assert r == pytest.approx([2.0, 0.707, 3.5355], rel=0.01) + + +def test_density_normalized(): + sample = numpy.random.normal(0, 1, (1000000, 2)) + queries = numpy.array([[0.0, 0.0], [-0.5, 0.7], [0.4, 1.7]]) + expected = numpy.exp(-(queries ** 2).sum(-1) / 2) / (2 * math.pi) + estimated = DTMDensity(k=150, normalize=True).fit(sample).transform(queries) + assert estimated == pytest.approx(expected, rel=0.4) + + +def test_density(): + distances = [[0, 1, 10], [2, 0, 30], [1, 3, 5]] + density = DTMDensity(k=2, metric="neighbors", dim=1).fit_transform(distances) + expected = numpy.array([2.0, 1.0, 0.5]) + assert density == pytest.approx(expected) + distances = [[0, 1], [2, 0], [1, 3]] + density = DTMDensity(metric="neighbors", dim=1).fit_transform(distances) + assert density == pytest.approx(expected) + density = DTMDensity(weights=[0.5, 0.5], metric="neighbors", dim=1).fit_transform(distances) + assert density == pytest.approx(expected) |