diff options
| author | Vincent Rouvreau <vincent.rouvreau@inria.fr> | 2022-01-21 08:24:16 +0100 |
|---|---|---|
| committer | Vincent Rouvreau <vincent.rouvreau@inria.fr> | 2022-01-21 08:24:16 +0100 |
| commit | 854cb0726336013cf9faede10ba61b0c6da938d3 (patch) | |
| tree | 2af6906f741904761064349f547811fb85b6000b /src/python | |
| parent | a6e7f96f7d2c391e4548309174cc05f5ae05d871 (diff) | |
| parent | de5aa9c891ef13c9fc2b2635bcd27ab873b0057b (diff) | |
Merge master
Diffstat (limited to 'src/python')
36 files changed, 1127 insertions, 204 deletions
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 73303a24..1314b444 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -14,13 +14,16 @@ function( add_GUDHI_PYTHON_lib THE_LIB ) endif(EXISTS ${THE_LIB}) endfunction( add_GUDHI_PYTHON_lib ) -function( add_GUDHI_PYTHON_lib_dir THE_LIB_DIR ) - # deals when it is not set - error on windows - if(EXISTS ${THE_LIB_DIR}) - set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${THE_LIB_DIR}', " PARENT_SCOPE) - else() - message("add_GUDHI_PYTHON_lib_dir - '${THE_LIB_DIR}' does not exist") - endif() +function( add_GUDHI_PYTHON_lib_dir) + # Argument may be a list (specifically on windows with release/debug paths) + foreach(THE_LIB_DIR IN LISTS ARGN) + # deals when it is not set - error on windows + if(EXISTS ${THE_LIB_DIR}) + set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${THE_LIB_DIR}', " PARENT_SCOPE) + else() + message("add_GUDHI_PYTHON_lib_dir - '${THE_LIB_DIR}' does not exist") + endif() + endforeach() endfunction( add_GUDHI_PYTHON_lib_dir ) # THE_TEST is the python test file name (without .py extension) containing tests functions @@ -41,13 +44,15 @@ function( add_gudhi_debug_info DEBUG_INFO ) endfunction( add_gudhi_debug_info ) if(PYTHONINTERP_FOUND) - if(PYBIND11_FOUND) + if(PYBIND11_FOUND AND CYTHON_FOUND) add_gudhi_debug_info("Pybind11 version ${PYBIND11_VERSION}") + # PyBind11 modules set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'bottleneck', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'hera', ") set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'clustering', ") - endif() - if(CYTHON_FOUND) + set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'datasets', ") + + # Cython modules set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'off_reader', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'simplex_tree', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'rips_complex', ") @@ -106,6 +111,16 @@ if(PYTHONINTERP_FOUND) if(TENSORFLOW_FOUND) add_gudhi_debug_info("TensorFlow version ${TENSORFLOW_VERSION}") endif() + if(SPHINX_FOUND) + add_gudhi_debug_info("Sphinx version ${SPHINX_VERSION}") + endif() + if(SPHINX_PARAMLINKS_FOUND) + add_gudhi_debug_info("Sphinx-paramlinks version ${SPHINX_PARAMLINKS_VERSION}") + endif() + if(PYTHON_DOCS_THEME_FOUND) + # Does not have a version number... + add_gudhi_debug_info("python_docs_theme found") + endif() set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DBOOST_RESULT_OF_USE_DECLTYPE', ") set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DBOOST_ALL_NO_LIB', ") @@ -151,6 +166,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera/wasserstein', ") set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera/bottleneck', ") if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'datasets/generators/_points', ") set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'bottleneck', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'nerve_gic', ") endif () @@ -163,6 +179,10 @@ if(PYTHONINTERP_FOUND) endif () if(CGAL_FOUND) + if(NOT CGAL_VERSION VERSION_LESS 5.3.0) + # CGAL_HEADER_ONLY has been dropped for CGAL >= 5.3. Only the header-only version is supported. + set(CGAL_HEADER_ONLY True) + endif(NOT CGAL_VERSION VERSION_LESS 5.3.0) # Add CGAL compilation args if(CGAL_HEADER_ONLY) add_gudhi_debug_info("CGAL header only version ${CGAL_VERSION}") @@ -170,7 +190,7 @@ if(PYTHONINTERP_FOUND) else(CGAL_HEADER_ONLY) add_gudhi_debug_info("CGAL version ${CGAL_VERSION}") add_GUDHI_PYTHON_lib("${CGAL_LIBRARY}") - add_GUDHI_PYTHON_lib_dir("${CGAL_LIBRARIES_DIR}") + add_GUDHI_PYTHON_lib_dir(${CGAL_LIBRARIES_DIR}) message("** Add CGAL ${CGAL_LIBRARIES_DIR}") # If CGAL is not header only, CGAL library may link with boost system, if(CMAKE_BUILD_TYPE MATCHES Debug) @@ -178,7 +198,7 @@ if(PYTHONINTERP_FOUND) else() add_GUDHI_PYTHON_lib("${Boost_SYSTEM_LIBRARY_RELEASE}") endif() - add_GUDHI_PYTHON_lib_dir("${Boost_LIBRARY_DIRS}") + add_GUDHI_PYTHON_lib_dir(${Boost_LIBRARY_DIRS}) message("** Add Boost ${Boost_LIBRARY_DIRS}") endif(CGAL_HEADER_ONLY) # GMP and GMPXX are not required, but if present, CGAL will link with them. @@ -190,13 +210,13 @@ if(PYTHONINTERP_FOUND) get_filename_component(GMP_LIBRARIES_DIR ${GMP_LIBRARIES} PATH) message("GMP_LIBRARIES_DIR from GMP_LIBRARIES set to ${GMP_LIBRARIES_DIR}") endif(NOT GMP_LIBRARIES_DIR) - add_GUDHI_PYTHON_lib_dir("${GMP_LIBRARIES_DIR}") + add_GUDHI_PYTHON_lib_dir(${GMP_LIBRARIES_DIR}) message("** Add gmp ${GMP_LIBRARIES_DIR}") if(GMPXX_FOUND) add_gudhi_debug_info("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}") set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_GMPXX', ") add_GUDHI_PYTHON_lib("${GMPXX_LIBRARIES}") - add_GUDHI_PYTHON_lib_dir("${GMPXX_LIBRARIES_DIR}") + add_GUDHI_PYTHON_lib_dir(${GMPXX_LIBRARIES_DIR}) message("** Add gmpxx ${GMPXX_LIBRARIES_DIR}") endif(GMPXX_FOUND) endif(GMP_FOUND) @@ -209,7 +229,7 @@ if(PYTHONINTERP_FOUND) get_filename_component(MPFR_LIBRARIES_DIR ${MPFR_LIBRARIES} PATH) message("MPFR_LIBRARIES_DIR from MPFR_LIBRARIES set to ${MPFR_LIBRARIES_DIR}") endif(NOT MPFR_LIBRARIES_DIR) - add_GUDHI_PYTHON_lib_dir("${MPFR_LIBRARIES_DIR}") + add_GUDHI_PYTHON_lib_dir(${MPFR_LIBRARIES_DIR}) message("** Add mpfr ${MPFR_LIBRARIES_DIR}") endif(MPFR_FOUND) endif(CGAL_FOUND) @@ -230,14 +250,14 @@ if(PYTHONINTERP_FOUND) if (TBB_FOUND AND WITH_GUDHI_USE_TBB) add_gudhi_debug_info("TBB version ${TBB_INTERFACE_VERSION} found and used") set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DGUDHI_USE_TBB', ") - if(CMAKE_BUILD_TYPE MATCHES Debug) + if((CMAKE_BUILD_TYPE MATCHES Debug) AND TBB_DEBUG_LIBRARY) add_GUDHI_PYTHON_lib("${TBB_DEBUG_LIBRARY}") add_GUDHI_PYTHON_lib("${TBB_MALLOC_DEBUG_LIBRARY}") else() add_GUDHI_PYTHON_lib("${TBB_RELEASE_LIBRARY}") add_GUDHI_PYTHON_lib("${TBB_MALLOC_RELEASE_LIBRARY}") endif() - add_GUDHI_PYTHON_lib_dir("${TBB_LIBRARY_DIRS}") + add_GUDHI_PYTHON_lib_dir(${TBB_LIBRARY_DIRS}) message("** Add tbb ${TBB_LIBRARY_DIRS}") set(GUDHI_PYTHON_INCLUDE_DIRS "${GUDHI_PYTHON_INCLUDE_DIRS}'${TBB_INCLUDE_DIRS}', ") endif() @@ -262,9 +282,12 @@ if(PYTHONINTERP_FOUND) file(COPY "gudhi/weighted_rips_complex.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") file(COPY "gudhi/dtm_rips_complex.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") file(COPY "gudhi/hera/__init__.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/hera") + file(COPY "gudhi/datasets" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi" FILES_MATCHING PATTERN "*.py") + # Some files for pip package file(COPY "introduction.rst" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/") + file(COPY "pyproject.toml" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/") add_custom_command( OUTPUT gudhi.so @@ -274,68 +297,74 @@ if(PYTHONINTERP_FOUND) add_custom_target(python ALL DEPENDS gudhi.so COMMENT "Do not forget to add ${CMAKE_CURRENT_BINARY_DIR}/ to your PYTHONPATH before using examples or tests") - install(CODE "execute_process(COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/setup.py install)") - - set(GUDHI_PYTHON_PATH_ENV "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}:$ENV{PYTHONPATH}") + # Path separator management for windows + if (WIN32) + set(GUDHI_PYTHON_PATH_ENV "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR};$ENV{PYTHONPATH}") + else(WIN32) + set(GUDHI_PYTHON_PATH_ENV "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}:$ENV{PYTHONPATH}") + endif(WIN32) # Documentation generation is available through sphinx - requires all modules # Make it first as sphinx test is by far the longest test which is nice when testing in parallel if(SPHINX_PATH) - if(MATPLOTLIB_FOUND) - if(NUMPY_FOUND) - if(SCIPY_FOUND) - if(SKLEARN_FOUND) - if(OT_FOUND) - if(PYBIND11_FOUND) - if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) - set (GUDHI_SPHINX_MESSAGE "Generating API documentation with Sphinx in ${CMAKE_CURRENT_BINARY_DIR}/sphinx/") - # User warning - Sphinx is a static pages generator, and configured to work fine with user_version - # Images and biblio warnings because not found on developper version - if (GUDHI_PYTHON_PATH STREQUAL "src/python") - set (GUDHI_SPHINX_MESSAGE "${GUDHI_SPHINX_MESSAGE} \n WARNING : Sphinx is configured for user version, you run it on developper version. Images and biblio will miss") - endif() - # sphinx target requires gudhi.so, because conf.py reads gudhi version from it - add_custom_target(sphinx - WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/doc - COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" - ${SPHINX_PATH} -b html ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/sphinx - DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/gudhi.so" - COMMENT "${GUDHI_SPHINX_MESSAGE}" VERBATIM) - - add_test(NAME sphinx_py_test - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" - ${SPHINX_PATH} -b doctest ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/doctest) - - # Set missing or not modules - set(GUDHI_MODULES ${GUDHI_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MODULES") - else(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) - message("++ Python documentation module will not be compiled because it requires a Eigen3 and CGAL version >= 4.11.0") + if(SPHINX_PARAMLINKS_FOUND) + if(PYTHON_DOCS_THEME_FOUND) + if(MATPLOTLIB_FOUND) + if(NUMPY_FOUND) + if(SCIPY_FOUND) + if(SKLEARN_FOUND) + if(OT_FOUND) + if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + set (GUDHI_SPHINX_MESSAGE "Generating API documentation with Sphinx in ${CMAKE_CURRENT_BINARY_DIR}/sphinx/") + # User warning - Sphinx is a static pages generator, and configured to work fine with user_version + # Images and biblio warnings because not found on developper version + if (GUDHI_PYTHON_PATH STREQUAL "src/python") + set (GUDHI_SPHINX_MESSAGE "${GUDHI_SPHINX_MESSAGE} \n WARNING : Sphinx is configured for user version, you run it on developper version. Images and biblio will miss") + endif() + # sphinx target requires gudhi.so, because conf.py reads gudhi version from it + add_custom_target(sphinx + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/doc + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" + ${SPHINX_PATH} -b html ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/sphinx + DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/gudhi.so" + COMMENT "${GUDHI_SPHINX_MESSAGE}" VERBATIM) + add_test(NAME sphinx_py_test + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" + ${SPHINX_PATH} -b doctest ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/doctest) + # Set missing or not modules + set(GUDHI_MODULES ${GUDHI_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MODULES") + else(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + message("++ Python documentation module will not be compiled because it requires a Eigen3 and CGAL version >= 4.11.0") + set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") + endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + else(OT_FOUND) + message("++ Python documentation module will not be compiled because POT was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) - else(PYBIND11_FOUND) - message("++ Python documentation module will not be compiled because pybind11 was not found") + endif(OT_FOUND) + else(SKLEARN_FOUND) + message("++ Python documentation module will not be compiled because scikit-learn was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(PYBIND11_FOUND) - else(OT_FOUND) - message("++ Python documentation module will not be compiled because POT was not found") + endif(SKLEARN_FOUND) + else(SCIPY_FOUND) + message("++ Python documentation module will not be compiled because scipy was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(OT_FOUND) - else(SKLEARN_FOUND) - message("++ Python documentation module will not be compiled because scikit-learn was not found") + endif(SCIPY_FOUND) + else(NUMPY_FOUND) + message("++ Python documentation module will not be compiled because numpy was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(SKLEARN_FOUND) - else(SCIPY_FOUND) - message("++ Python documentation module will not be compiled because scipy was not found") + endif(NUMPY_FOUND) + else(MATPLOTLIB_FOUND) + message("++ Python documentation module will not be compiled because matplotlib was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(SCIPY_FOUND) - else(NUMPY_FOUND) - message("++ Python documentation module will not be compiled because numpy was not found") + endif(MATPLOTLIB_FOUND) + else(PYTHON_DOCS_THEME_FOUND) + message("++ Python documentation module will not be compiled because python-docs-theme was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(NUMPY_FOUND) - else(MATPLOTLIB_FOUND) - message("++ Python documentation module will not be compiled because matplotlib was not found") + endif(PYTHON_DOCS_THEME_FOUND) + else(SPHINX_PARAMLINKS_FOUND) + message("++ Python documentation module will not be compiled because sphinxcontrib-paramlinks was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(MATPLOTLIB_FOUND) + endif(SPHINX_PARAMLINKS_FOUND) else(SPHINX_PATH) message("++ Python documentation module will not be compiled because sphinx and sphinxcontrib-bibtex were not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") @@ -383,9 +412,7 @@ if(PYTHONINTERP_FOUND) COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/bottleneck_basic_example.py") - if (PYBIND11_FOUND) - add_gudhi_py_test(test_bottleneck_distance) - endif() + add_gudhi_py_test(test_bottleneck_distance) # Cover complex file(COPY ${CMAKE_SOURCE_DIR}/data/points/human.off DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) @@ -425,6 +452,10 @@ if(PYTHONINTERP_FOUND) WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/alpha_complex_from_points_example.py") + add_test(NAME alpha_complex_from_generated_points_on_sphere_example_py_test + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" + ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/alpha_complex_from_generated_points_on_sphere_example.py") add_test(NAME alpha_complex_diagram_persistence_from_off_file_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" @@ -437,6 +468,9 @@ if(PYTHONINTERP_FOUND) # Euclidean witness add_gudhi_py_test(test_euclidean_witness_complex) + # Datasets generators + add_gudhi_py_test(test_datasets_generators) # TODO separate full python datasets generators in another test file independant from CGAL ? + endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) # Cubical @@ -459,7 +493,7 @@ if(PYTHONINTERP_FOUND) WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py" - --no-diagram -f ${CMAKE_SOURCE_DIR}/data/distance_matrix/lower_triangular_distance_matrix.csv -e 12.0 -d 3) + --no-diagram -f ${CMAKE_SOURCE_DIR}/data/distance_matrix/lower_triangular_distance_matrix.csv -s , -e 12.0 -d 3) add_test(NAME rips_complex_diagram_persistence_from_off_file_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} @@ -495,14 +529,14 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_reader_utils) # Wasserstein - if(OT_FOUND AND PYBIND11_FOUND) + if(OT_FOUND) # EagerPy dependency because of enable_autodiff=True if(EAGERPY_FOUND) add_gudhi_py_test(test_wasserstein_distance) endif() + add_gudhi_py_test(test_wasserstein_barycenter) - endif() - if(OT_FOUND) + if(TORCH_FOUND AND TENSORFLOW_FOUND AND EAGERPY_FOUND) add_gudhi_py_test(test_wasserstein_with_tensors) endif() @@ -523,7 +557,7 @@ if(PYTHONINTERP_FOUND) endif() # Tomato - if(SCIPY_FOUND AND SKLEARN_FOUND AND PYBIND11_FOUND) + if(SCIPY_FOUND AND SKLEARN_FOUND) add_gudhi_py_test(test_tomato) endif() @@ -540,10 +574,10 @@ if(PYTHONINTERP_FOUND) # Set missing or not modules set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES") - else(CYTHON_FOUND) - message("++ Python module will not be compiled because cython was not found") + else(PYBIND11_FOUND AND CYTHON_FOUND) + message("++ Python module will not be compiled because cython and/or pybind11 was/were not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(CYTHON_FOUND) + endif(PYBIND11_FOUND AND CYTHON_FOUND) else(PYTHONINTERP_FOUND) message("++ Python module will not be compiled because no Python interpreter was found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python" CACHE INTERNAL "GUDHI_MISSING_MODULES") diff --git a/src/python/doc/_templates/layout.html b/src/python/doc/_templates/layout.html index cd40a51b..e074b6c7 100644 --- a/src/python/doc/_templates/layout.html +++ b/src/python/doc/_templates/layout.html @@ -194,6 +194,7 @@ <li><a href="/relatedprojects/">Related projects</a></li> <li><a href="/theyaretalkingaboutus/">They are talking about us</a></li> <li><a href="/inaction/">GUDHI in action</a></li> + <li><a href="/etymology/">Etymology</a></li> </ul> </li> <li class="divider"></li> diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst index fffcb3db..96e267ef 100644 --- a/src/python/doc/alpha_complex_user.rst +++ b/src/python/doc/alpha_complex_user.rst @@ -163,7 +163,10 @@ As the squared radii computed by CGAL are an approximation, it might happen that :math:`\alpha^2` values do not quite define a proper filtration (i.e. non-decreasing with respect to inclusion). We fix that up by calling :func:`~gudhi.SimplexTree.make_filtration_non_decreasing` (cf. -`C++ version <http://gudhi.gforge.inria.fr/doc/latest/index.html>`_). +`C++ version <https://gudhi.inria.fr/doc/latest/class_gudhi_1_1_simplex__tree.html>`_). + +.. note:: + This is not the case in `exact` version, this is the reason why it is not called in this case. Prune above given filtration value ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/src/python/doc/conf.py b/src/python/doc/conf.py index b06baf9c..e69e2751 100755 --- a/src/python/doc/conf.py +++ b/src/python/doc/conf.py @@ -120,15 +120,12 @@ pygments_style = 'sphinx' # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. -html_theme = 'classic' +html_theme = 'python_docs_theme' # Theme options are theme-specific and customize the look and feel of a theme # further. For a list of options available for each theme, see the # documentation. html_theme_options = { - "sidebarbgcolor": "#A1ADCD", - "sidebartextcolor": "black", - "sidebarlinkcolor": "#334D5C", "body_max_width": "100%", } diff --git a/src/python/doc/datasets_generators.inc b/src/python/doc/datasets_generators.inc new file mode 100644 index 00000000..8d169275 --- /dev/null +++ b/src/python/doc/datasets_generators.inc @@ -0,0 +1,14 @@ +.. table:: + :widths: 30 40 30 + + +-----------------------------------+--------------------------------------------+--------------------------------------------------------------------------------------+ + | .. figure:: | Datasets generators (points). | :Authors: Hind Montassif | + | img/sphere_3d.png | | | + | | | :Since: GUDHI 3.5.0 | + | | | | + | | | :License: MIT (`LGPL v3 </licensing/>`_) | + | | | | + | | | :Requires: `CGAL <installation.html#cgal>`_ | + +-----------------------------------+--------------------------------------------+--------------------------------------------------------------------------------------+ + | * :doc:`datasets_generators` | + +-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/doc/datasets_generators.rst b/src/python/doc/datasets_generators.rst new file mode 100644 index 00000000..260c3882 --- /dev/null +++ b/src/python/doc/datasets_generators.rst @@ -0,0 +1,105 @@ + +:orphan: + +.. To get rid of WARNING: document isn't included in any toctree + +=========================== +Datasets generators manual +=========================== + +We provide the generation of different customizable datasets to use as inputs for Gudhi complexes and data structures. + + +Points generators +------------------ + +The module **points** enables the generation of random points on a sphere, random points on a torus and as a grid. + +Points on sphere +^^^^^^^^^^^^^^^^ + +The function **sphere** enables the generation of random i.i.d. points uniformly on a (d-1)-sphere in :math:`R^d`. +The user should provide the number of points to be generated on the sphere :code:`n_samples` and the ambient dimension :code:`ambient_dim`. +The :code:`radius` of sphere is optional and is equal to **1** by default. +Only random points generation is currently available. + +The generated points are given as an array of shape :math:`(n\_samples, ambient\_dim)`. + +Example +""""""" + +.. code-block:: python + + from gudhi.datasets.generators import points + from gudhi import AlphaComplex + + # Generate 50 points on a sphere in R^2 + gen_points = points.sphere(n_samples = 50, ambient_dim = 2, radius = 1, sample = "random") + + # Create an alpha complex from the generated points + alpha_complex = AlphaComplex(points = gen_points) + +.. autofunction:: gudhi.datasets.generators.points.sphere + +Points on a flat torus +^^^^^^^^^^^^^^^^^^^^^^ + +You can also generate points on a torus. + +Two functions are available and give the same output: the first one depends on **CGAL** and the second does not and consists of full python code. + +On another hand, two sample types are provided: you can either generate i.i.d. points on a d-torus in :math:`R^{2d}` *randomly* or on a *grid*. + +First function: **ctorus** +""""""""""""""""""""""""""" + +The user should provide the number of points to be generated on the torus :code:`n_samples`, and the dimension :code:`dim` of the torus on which points would be generated in :math:`R^{2dim}`. +The :code:`sample` argument is optional and is set to **'random'** by default. +In this case, the returned generated points would be an array of shape :math:`(n\_samples, 2*dim)`. +Otherwise, if set to **'grid'**, the points are generated on a grid and would be given as an array of shape: + +.. math:: + + ( ⌊n\_samples^{1 \over {dim}}⌋^{dim}, 2*dim ) + +**Note 1:** The output array first shape is rounded down to the closest perfect :math:`dim^{th}` power. + +**Note 2:** This version is recommended when the user wishes to use **'grid'** as sample type, or **'random'** with a relatively small number of samples (~ less than 150). + +Example +""""""" +.. code-block:: python + + from gudhi.datasets.generators import points + + # Generate 50 points randomly on a torus in R^6 + gen_points = points.ctorus(n_samples = 50, dim = 3) + + # Generate 27 points on a torus as a grid in R^6 + gen_points = points.ctorus(n_samples = 50, dim = 3, sample = 'grid') + +.. autofunction:: gudhi.datasets.generators.points.ctorus + +Second function: **torus** +""""""""""""""""""""""""""" + +The user should provide the number of points to be generated on the torus :code:`n_samples` and the dimension :code:`dim` of the torus on which points would be generated in :math:`R^{2dim}`. +The :code:`sample` argument is optional and is set to **'random'** by default. +The other allowed value of sample type is **'grid'**. + +**Note:** This version is recommended when the user wishes to use **'random'** as sample type with a great number of samples and a low dimension. + +Example +""""""" +.. code-block:: python + + from gudhi.datasets.generators import points + + # Generate 50 points randomly on a torus in R^6 + gen_points = points.torus(n_samples = 50, dim = 3) + + # Generate 27 points on a torus as a grid in R^6 + gen_points = points.torus(n_samples = 50, dim = 3, sample = 'grid') + + +.. autofunction:: gudhi.datasets.generators.points.torus diff --git a/src/python/doc/examples.rst b/src/python/doc/examples.rst index 76e5d4c7..1442f185 100644 --- a/src/python/doc/examples.rst +++ b/src/python/doc/examples.rst @@ -8,6 +8,7 @@ Examples .. only:: builder_html * :download:`alpha_complex_diagram_persistence_from_off_file_example.py <../example/alpha_complex_diagram_persistence_from_off_file_example.py>` + * :download:`alpha_complex_from_generated_points_on_sphere_example.py <../example/alpha_complex_from_generated_points_on_sphere_example.py>` * :download:`alpha_complex_from_points_example.py <../example/alpha_complex_from_points_example.py>` * :download:`alpha_rips_persistence_bottleneck_distance.py <../example/alpha_rips_persistence_bottleneck_distance.py>` * :download:`bottleneck_basic_example.py <../example/bottleneck_basic_example.py>` diff --git a/src/python/doc/img/sphere_3d.png b/src/python/doc/img/sphere_3d.png Binary files differnew file mode 100644 index 00000000..70f3184f --- /dev/null +++ b/src/python/doc/img/sphere_3d.png diff --git a/src/python/doc/index.rst b/src/python/doc/index.rst index 040e57a4..2d7921ae 100644 --- a/src/python/doc/index.rst +++ b/src/python/doc/index.rst @@ -91,3 +91,8 @@ Clustering ********** .. include:: clustering.inc + +Datasets generators +******************* + +.. include:: datasets_generators.inc diff --git a/src/python/doc/installation.rst b/src/python/doc/installation.rst index 66efe45a..35c344e3 100644 --- a/src/python/doc/installation.rst +++ b/src/python/doc/installation.rst @@ -41,7 +41,7 @@ there. The library uses c++14 and requires `Boost <https://www.boost.org/>`_ :math:`\geq` 1.56.0, `CMake <https://www.cmake.org/>`_ :math:`\geq` 3.5 to generate makefiles, -`NumPy <http://numpy.org>`_, `Cython <https://www.cython.org/>`_ and +`NumPy <http://numpy.org>`_ :math:`\geq` 1.15.0, `Cython <https://www.cython.org/>`_ and `pybind11 <https://github.com/pybind/pybind11>`_ to compile the GUDHI Python module. It is a multi-platform library and compiles on Linux, Mac OSX and Visual @@ -99,20 +99,14 @@ Or install it definitely in your Python packages folder: .. code-block:: bash cd /path-to-gudhi/build/python - # May require sudo or administrator privileges - make install + python setup.py install # add --user to the command if you do not have the permission + # Or 'pip install .' .. note:: - :code:`make install` is only a - `CMake custom targets <https://cmake.org/cmake/help/latest/command/add_custom_target.html>`_ - to shortcut :code:`python setup.py install` command. It does not take into account :code:`CMAKE_INSTALL_PREFIX`. - But one can use :code:`python setup.py install ...` specific options in the python directory: - -.. code-block:: bash - - python setup.py install --prefix /home/gudhi # Install in /home/gudhi directory + But one can use + `alternate location installation <https://docs.python.org/3/install/#alternate-installation>`_. Test suites =========== @@ -200,8 +194,10 @@ A complete configuration would be : Documentation ============= -To build the documentation, `sphinx-doc <http://www.sphinx-doc.org>`_ and -`sphinxcontrib-bibtex <https://sphinxcontrib-bibtex.readthedocs.io>`_ are +To build the documentation, `sphinx-doc <http://www.sphinx-doc.org>`_, +`sphinxcontrib-bibtex <https://sphinxcontrib-bibtex.readthedocs.io>`_, +`sphinxcontrib-paramlinks <https://github.com/sqlalchemyorg/sphinx-paramlinks>`_ and +`python-docs-theme <https://github.com/python/python-docs-theme>`_ are required. As the documentation is auto-tested, `CGAL`_, `Eigen`_, `Matplotlib`_, `NumPy`_, `POT`_, `Scikit-learn`_ and `SciPy`_ are also mandatory to build the documentation. @@ -363,7 +359,7 @@ Python Optimal Transport ------------------------ The :doc:`Wasserstein distance </wasserstein_distance_user>` -module requires `POT <https://pot.readthedocs.io/>`_, a library that provides +module requires `POT <https://pythonot.github.io/>`_, a library that provides several solvers for optimization problems related to Optimal Transport. PyTorch @@ -402,8 +398,9 @@ TensorFlow Bug reports and contributions ***************************** -Please help us improving the quality of the GUDHI library. You may report bugs or suggestions to: - - Contact: gudhi-users@lists.gforge.inria.fr +Please help us improving the quality of the GUDHI library. +You may `report bugs <https://github.com/GUDHI/gudhi-devel/issues>`_ or +`contact us <https://gudhi.inria.fr/contact/>`_ for any suggestions. -GUDHI is open to external contributions. If you want to join our development team, please contact us. +GUDHI is open to external contributions. If you want to join our development team, please take some time to read our +`contributing guide <https://github.com/GUDHI/gudhi-devel/blob/master/.github/CONTRIBUTING.md>`_. diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index 9ffc2759..76eb1469 100644 --- a/src/python/doc/wasserstein_distance_user.rst +++ b/src/python/doc/wasserstein_distance_user.rst @@ -44,7 +44,7 @@ Basic example ************* This example computes the 1-Wasserstein distance from 2 persistence diagrams with Euclidean ground metric. -Note that persistence diagrams must be submitted as (n x 2) numpy arrays and must not contain inf values. +Note that persistence diagrams must be submitted as (n x 2) numpy arrays. .. testcode:: @@ -67,14 +67,16 @@ We can also have access to the optimal matching by letting `matching=True`. It is encoded as a list of indices (i,j), meaning that the i-th point in X is mapped to the j-th point in Y. An index of -1 represents the diagonal. +It handles essential parts (points with infinite coordinates). However if the cardinalities of the essential parts differ, +any matching has a cost +inf and thus can be considered to be optimal. In such a case, the function returns `(np.inf, None)`. .. testcode:: import gudhi.wasserstein import numpy as np - dgm1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]]) - dgm2 = np.array([[2.8, 4.45], [5, 6], [9.5, 14.1]]) + dgm1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974], [3, np.inf]]) + dgm2 = np.array([[2.8, 4.45], [5, 6], [9.5, 14.1], [4, np.inf]]) cost, matchings = gudhi.wasserstein.wasserstein_distance(dgm1, dgm2, matching=True, order=1, internal_p=2) message_cost = "Wasserstein distance value = %.2f" %cost @@ -90,16 +92,31 @@ An index of -1 represents the diagonal. for j in dgm2_to_diagonal: print("point %s in dgm2 is matched to the diagonal" %j) -The output is: + # An example where essential part cardinalities differ + dgm3 = np.array([[1, 2], [0, np.inf]]) + dgm4 = np.array([[1, 2], [0, np.inf], [1, np.inf]]) + cost, matchings = gudhi.wasserstein.wasserstein_distance(dgm3, dgm4, matching=True, order=1, internal_p=2) + print("\nSecond example:") + print("cost:", cost) + print("matchings:", matchings) + + +The output is: .. testoutput:: - Wasserstein distance value = 2.15 + Wasserstein distance value = 3.15 point 0 in dgm1 is matched to point 0 in dgm2 point 1 in dgm1 is matched to point 2 in dgm2 + point 3 in dgm1 is matched to point 3 in dgm2 point 2 in dgm1 is matched to the diagonal point 1 in dgm2 is matched to the diagonal + Second example: + cost: inf + matchings: None + + Barycenters ----------- @@ -181,4 +198,4 @@ Tutorial This `notebook <https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-Barycenters-of-persistence-diagrams.ipynb>`_ -presents the concept of barycenter, or Fréchet mean, of a family of persistence diagrams.
\ No newline at end of file +presents the concept of barycenter, or Fréchet mean, of a family of persistence diagrams. diff --git a/src/python/example/alpha_complex_from_generated_points_on_sphere_example.py b/src/python/example/alpha_complex_from_generated_points_on_sphere_example.py new file mode 100644 index 00000000..3558077e --- /dev/null +++ b/src/python/example/alpha_complex_from_generated_points_on_sphere_example.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python + +from gudhi.datasets.generators import _points +from gudhi import AlphaComplex + + +""" 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 +""" + +__author__ = "Hind Montassif" +__copyright__ = "Copyright (C) 2021 Inria" +__license__ = "MIT" + +print("#####################################################################") +print("AlphaComplex creation from generated points on sphere") + + +gen_points = _points.sphere(n_samples = 50, ambient_dim = 2, radius = 1, sample = "random") + +# Create an alpha complex +alpha_complex = AlphaComplex(points = gen_points) +simplex_tree = alpha_complex.create_simplex_tree() + +result_str = 'Alpha complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ + repr(simplex_tree.num_simplices()) + ' simplices - ' + \ + repr(simplex_tree.num_vertices()) + ' vertices.' +print(result_str) + diff --git a/src/python/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py b/src/python/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py index 236d085d..8a9cc857 100755 --- a/src/python/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py +++ b/src/python/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py @@ -21,11 +21,12 @@ parser = argparse.ArgumentParser( description="RipsComplex creation from " "a distance matrix read in a csv file.", epilog="Example: " "example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py " - "-f ../data/distance_matrix/lower_triangular_distance_matrix.csv -e 12.0 -d 3" + "-f ../data/distance_matrix/lower_triangular_distance_matrix.csv -s , -e 12.0 -d 3" "- Constructs a Rips complex with the " "distance matrix from the given csv file.", ) parser.add_argument("-f", "--file", type=str, required=True) +parser.add_argument("-s", "--separator", type=str, required=True) parser.add_argument("-e", "--max_edge_length", type=float, default=0.5) parser.add_argument("-d", "--max_dimension", type=int, default=1) parser.add_argument("-b", "--band", type=float, default=0.0) @@ -44,7 +45,7 @@ print("RipsComplex creation from distance matrix read in a csv file") message = "RipsComplex with max_edge_length=" + repr(args.max_edge_length) print(message) -distance_matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file=args.file) +distance_matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file=args.file, separator=args.separator) rips_complex = gudhi.RipsComplex( distance_matrix=distance_matrix, max_edge_length=args.max_edge_length ) diff --git a/src/python/gudhi/clustering/tomato.py b/src/python/gudhi/clustering/tomato.py index fbba3cc8..d0e9995c 100644 --- a/src/python/gudhi/clustering/tomato.py +++ b/src/python/gudhi/clustering/tomato.py @@ -271,7 +271,7 @@ class Tomato: l = self.max_weight_per_cc_.min() r = self.max_weight_per_cc_.max() if self.diagram_.size > 0: - plt.plot(self.diagram_[:, 0], self.diagram_[:, 1], "ro") + plt.plot(self.diagram_[:, 0], self.diagram_[:, 1], "o", color="red") l = min(l, self.diagram_[:, 1].min()) r = max(r, self.diagram_[:, 0].max()) if l == r: @@ -283,7 +283,7 @@ class Tomato: l, r = -1.0, 1.0 plt.plot([l, r], [l, r]) plt.plot( - self.max_weight_per_cc_, numpy.full(self.max_weight_per_cc_.shape, 1.1 * l - 0.1 * r), "ro", color="green" + self.max_weight_per_cc_, numpy.full(self.max_weight_per_cc_.shape, 1.1 * l - 0.1 * r), "o", color="green" ) plt.show() 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/__init__.py b/src/python/gudhi/datasets/__init__.py new file mode 100644 index 00000000..e69de29b --- /dev/null +++ b/src/python/gudhi/datasets/__init__.py diff --git a/src/python/gudhi/datasets/generators/__init__.py b/src/python/gudhi/datasets/generators/__init__.py new file mode 100644 index 00000000..e69de29b --- /dev/null +++ b/src/python/gudhi/datasets/generators/__init__.py 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.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 994be3b6..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 @@ -111,7 +112,7 @@ class KNearestNeighbors: nargs = { k: v for k, v in self.params.items() if k in {"p", "n_jobs", "metric_params", "algorithm", "leaf_size"} } - self.nn = NearestNeighbors(self.k, metric=self.metric, **nargs) + self.nn = NearestNeighbors(n_neighbors=self.k, metric=self.metric, **nargs) self.nn.fit(X) if self.params["implementation"] == "hnsw": @@ -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 3df614dd..70311ead 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -44,7 +44,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_interface_full_featured "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>": - Simplex_tree() nogil + Simplex_tree_interface_full_featured() nogil double simplex_filtration(vector[int] simplex) nogil void assign_simplex_filtration(vector[int] simplex, double filtration) nogil void initialize_filtration() nogil @@ -79,7 +79,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 67428401..6393343f 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -9,9 +9,8 @@ from cython.operator import dereference, preincrement from libc.stdint cimport intptr_t -import numpy -from numpy import array as np_array -cimport simplex_tree +import numpy as np +cimport gudhi.simplex_tree __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" @@ -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): diff --git a/src/python/gudhi/wasserstein/wasserstein.py b/src/python/gudhi/wasserstein/wasserstein.py index a9d1cdff..dc18806e 100644 --- a/src/python/gudhi/wasserstein/wasserstein.py +++ b/src/python/gudhi/wasserstein/wasserstein.py @@ -9,6 +9,7 @@ import numpy as np import scipy.spatial.distance as sc +import warnings try: import ot @@ -70,6 +71,7 @@ def _perstot_autodiff(X, order, internal_p): ''' return _dist_to_diag(X, internal_p).norms.lp(order) + def _perstot(X, order, internal_p, enable_autodiff): ''' :param X: (n x 2) numpy.array (points of a given diagram). @@ -79,6 +81,9 @@ def _perstot(X, order, internal_p, enable_autodiff): transparent to automatic differentiation. :type enable_autodiff: bool :returns: float, the total persistence of the diagram (that is, its distance to the empty diagram). + + .. note:: + Can be +inf if the diagram has an essential part (points with infinite coordinates). ''' if enable_autodiff: import eagerpy as ep @@ -88,32 +93,163 @@ def _perstot(X, order, internal_p, enable_autodiff): return np.linalg.norm(_dist_to_diag(X, internal_p), ord=order) -def wasserstein_distance(X, Y, matching=False, order=1., internal_p=np.inf, enable_autodiff=False): +def _get_essential_parts(a): ''' - :param X: (n x 2) numpy.array encoding the (finite points of the) first diagram. Must not contain essential points - (i.e. with infinite coordinate). - :param Y: (m x 2) numpy.array encoding the second diagram. - :param matching: if True, computes and returns the optimal matching between X and Y, encoded as - a (n x 2) np.array [...[i,j]...], meaning the i-th point in X is matched to - the j-th point in Y, with the convention (-1) represents the diagonal. - :param order: exponent for Wasserstein; Default value is 1. - :param internal_p: Ground metric on the (upper-half) plane (i.e. norm L^p in R^2); - Default value is `np.inf`. - :param enable_autodiff: If X and Y are torch.tensor or tensorflow.Tensor, make the computation + :param a: (n x 2) numpy.array (point of a diagram) + :returns: five lists of indices (between 0 and len(a)) accounting for the five types of points with infinite + coordinates that can occur in a diagram, namely: + type0 : (-inf, finite) + type1 : (finite, +inf) + type2 : (-inf, +inf) + type3 : (-inf, -inf) + type4 : (+inf, +inf) + .. note:: + For instance, a[_get_essential_parts(a)[0]] returns the points in a of coordinates (-inf, x) for some finite x. + Note also that points with (+inf, -inf) are not handled (points (x,y) in dgm satisfy by assumption (y >= x)). + + Finally, we consider that points with coordinates (-inf,-inf) and (+inf, +inf) belong to the diagonal. + ''' + if len(a): + first_coord_finite = np.isfinite(a[:,0]) + second_coord_finite = np.isfinite(a[:,1]) + first_coord_infinite_positive = (a[:,0] == np.inf) + second_coord_infinite_positive = (a[:,1] == np.inf) + first_coord_infinite_negative = (a[:,0] == -np.inf) + second_coord_infinite_negative = (a[:,1] == -np.inf) + + ess_first_type = np.where(second_coord_finite & first_coord_infinite_negative)[0] # coord (-inf, x) + ess_second_type = np.where(first_coord_finite & second_coord_infinite_positive)[0] # coord (x, +inf) + ess_third_type = np.where(first_coord_infinite_negative & second_coord_infinite_positive)[0] # coord (-inf, +inf) + + ess_fourth_type = np.where(first_coord_infinite_negative & second_coord_infinite_negative)[0] # coord (-inf, -inf) + ess_fifth_type = np.where(first_coord_infinite_positive & second_coord_infinite_positive)[0] # coord (+inf, +inf) + return ess_first_type, ess_second_type, ess_third_type, ess_fourth_type, ess_fifth_type + else: + return [], [], [], [], [] + + +def _cost_and_match_essential_parts(X, Y, idX, idY, order, axis): + ''' + :param X: (n x 2) numpy.array (dgm points) + :param Y: (n x 2) numpy.array (dgm points) + :param idX: indices to consider for this one dimensional OT problem (in X) + :param idY: indices to consider for this one dimensional OT problem (in Y) + :param order: exponent for Wasserstein distance computation + :param axis: must be 0 or 1, correspond to the coordinate which is finite. + :returns: cost (float) and match for points with *one* infinite coordinate. + + .. note:: + Assume idX, idY come when calling _handle_essential_parts, thus have same length. + ''' + u = X[idX, axis] + v = Y[idY, axis] + + cost = np.sum(np.abs(np.sort(u) - np.sort(v))**(order)) # OT cost in 1D + + sortidX = idX[np.argsort(u)] + sortidY = idY[np.argsort(v)] + # We return [i,j] sorted per value + match = list(zip(sortidX, sortidY)) + + return cost, match + + +def _handle_essential_parts(X, Y, order): + ''' + :param X: (n x 2) numpy array, first diagram. + :param Y: (n x 2) numpy array, second diagram. + :order: Wasserstein order for cost computation. + :returns: cost and matching due to essential parts. If cost is +inf, matching will be set to None. + ''' + ess_parts_X = _get_essential_parts(X) + ess_parts_Y = _get_essential_parts(Y) + + # Treats the case of infinite cost (cardinalities of essential parts differ). + for u, v in list(zip(ess_parts_X, ess_parts_Y))[:3]: # ignore types 4 and 5 as they belong to the diagonal + if len(u) != len(v): + return np.inf, None + + # Now we know each essential part has the same number of points in both diagrams. + # Handle type 0 and type 1 essential parts (those with one finite coordinates) + c1, m1 = _cost_and_match_essential_parts(X, Y, ess_parts_X[0], ess_parts_Y[0], axis=1, order=order) + c2, m2 = _cost_and_match_essential_parts(X, Y, ess_parts_X[1], ess_parts_Y[1], axis=0, order=order) + + c = c1 + c2 + m = m1 + m2 + + # Handle type3 (coordinates (-inf,+inf), so we just align points) + m += list(zip(ess_parts_X[2], ess_parts_Y[2])) + + # Handle type 4 and 5, considered as belonging to the diagonal so matched to (-1) with cost 0. + for z in ess_parts_X[3:]: + m += [(u, -1) for u in z] # points in X are matched to -1 + for z in ess_parts_Y[3:]: + m += [(-1, v) for v in z] # -1 is match to points in Y + + return c, np.array(m) + + +def _finite_part(X): + ''' + :param X: (n x 2) numpy array encoding a persistence diagram. + :returns: The finite part of a diagram `X` (points with finite coordinates). + ''' + return X[np.where(np.isfinite(X[:,0]) & np.isfinite(X[:,1]))] + + +def _warn_infty(matching): + ''' + Handle essential parts with different cardinalities. Warn the user about cost being infinite and (if + `matching=True`) about the returned matching being `None`. + ''' + if matching: + warnings.warn('Cardinality of essential parts differs. Distance (cost) is +inf, and the returned matching is None.') + return np.inf, None + else: + warnings.warn('Cardinality of essential parts differs. Distance (cost) is +inf.') + return np.inf + + +def wasserstein_distance(X, Y, matching=False, order=1., internal_p=np.inf, enable_autodiff=False, + keep_essential_parts=True): + ''' + Compute the Wasserstein distance between persistence diagram using Python Optimal Transport backend. + Diagrams can contain points with infinity coordinates (essential parts). + Points with (-inf,-inf) and (+inf,+inf) coordinates are considered as belonging to the diagonal. + If the distance between two diagrams is +inf (which happens if the cardinalities of essential + parts differ) and optimal matching is required, it will be set to ``None``. + + :param X: The first diagram. + :type X: n x 2 numpy.array + :param Y: The second diagram. + :type Y: m x 2 numpy.array + :param matching: if ``True``, computes and returns the optimal matching between X and Y, encoded as + a (n x 2) np.array [...[i,j]...], meaning the i-th point in X is matched to + the j-th point in Y, with the convention that (-1) represents the diagonal. + :param order: Wasserstein exponent q (1 <= q < infinity). + :type order: float + :param internal_p: Ground metric on the (upper-half) plane (i.e. norm L^p in R^2). + :type internal_p: float + :param enable_autodiff: If X and Y are ``torch.tensor`` or ``tensorflow.Tensor``, make the computation transparent to automatic differentiation. This requires the package EagerPy and is currently incompatible - with `matching=True`. + with ``matching=True`` and with ``keep_essential_parts=True``. - .. note:: This considers the function defined on the coordinates of the off-diagonal points of X and Y + .. note:: This considers the function defined on the coordinates of the off-diagonal finite points of X and Y and lets the various frameworks compute its gradient. It never pulls new points from the diagonal. :type enable_autodiff: bool - :returns: the Wasserstein distance of order q (1 <= q < infinity) between persistence diagrams with + :param keep_essential_parts: If ``False``, only considers the finite points in the diagrams. + Otherwise, include essential parts in cost and matching computation. + :type keep_essential_parts: bool + :returns: The Wasserstein distance of order q (1 <= q < infinity) between persistence diagrams with respect to the internal_p-norm as ground metric. If matching is set to True, also returns the optimal matching between X and Y. + If cost is +inf, any matching is optimal and thus it returns `None` instead. ''' + + # First step: handle empty diagrams n = len(X) m = len(Y) - # handle empty diagrams if n == 0: if m == 0: if not matching: @@ -122,16 +258,45 @@ def wasserstein_distance(X, Y, matching=False, order=1., internal_p=np.inf, enab else: return 0., np.array([]) else: - if not matching: - return _perstot(Y, order, internal_p, enable_autodiff) + cost = _perstot(Y, order, internal_p, enable_autodiff) + if cost == np.inf: + return _warn_infty(matching) else: - return _perstot(Y, order, internal_p, enable_autodiff), np.array([[-1, j] for j in range(m)]) + if not matching: + return cost + else: + return cost, np.array([[-1, j] for j in range(m)]) elif m == 0: - if not matching: - return _perstot(X, order, internal_p, enable_autodiff) + cost = _perstot(X, order, internal_p, enable_autodiff) + if cost == np.inf: + return _warn_infty(matching) else: - return _perstot(X, order, internal_p, enable_autodiff), np.array([[i, -1] for i in range(n)]) + if not matching: + return cost + else: + return cost, np.array([[i, -1] for i in range(n)]) + + # Check essential part and enable autodiff together + if enable_autodiff and keep_essential_parts: + warnings.warn('''enable_autodiff=True and keep_essential_parts=True are incompatible together. + keep_essential_parts is set to False: only points with finite coordinates are considered + in the following. + ''') + keep_essential_parts = False + + # Second step: handle essential parts if needed. + if keep_essential_parts: + essential_cost, essential_matching = _handle_essential_parts(X, Y, order=order) + if (essential_cost == np.inf): + return _warn_infty(matching) # Tells the user that cost is infty and matching (if True) is None. + # avoid computing transport cost between the finite parts if essential parts + # cardinalities do not match (saves time) + else: + essential_cost = 0 + essential_matching = None + + # Now the standard pipeline for finite parts if enable_autodiff: import eagerpy as ep @@ -139,6 +304,12 @@ def wasserstein_distance(X, Y, matching=False, order=1., internal_p=np.inf, enab Y_orig = ep.astensor(Y) X = X_orig.numpy() Y = Y_orig.numpy() + + # Extract finite points of the diagrams. + X, Y = _finite_part(X), _finite_part(Y) + n = len(X) + m = len(Y) + M = _build_dist_matrix(X, Y, order=order, internal_p=internal_p) a = np.ones(n+1) # weight vector of the input diagram. Uniform here. a[-1] = m @@ -154,7 +325,10 @@ def wasserstein_distance(X, Y, matching=False, order=1., internal_p=np.inf, enab # Now we turn to -1 points encoding the diagonal match[:,0][match[:,0] >= n] = -1 match[:,1][match[:,1] >= m] = -1 - return ot_cost ** (1./order) , match + # Finally incorporate the essential part matching + if essential_matching is not None: + match = np.concatenate([match, essential_matching]) if essential_matching.size else match + return (ot_cost + essential_cost) ** (1./order) , match if enable_autodiff: P = ot.emd(a=a, b=b, M=M, numItermax=2000000) @@ -173,9 +347,9 @@ def wasserstein_distance(X, Y, matching=False, order=1., internal_p=np.inf, enab return ep.concatenate(dists).norms.lp(order).raw # We can also concatenate the 3 vectors to compute just one norm. - # Comptuation of the otcost using the ot.emd2 library. + # Comptuation of the ot cost using the ot.emd2 library. # Note: it is the Wasserstein distance to the power q. # The default numItermax=100000 is not sufficient for some examples with 5000 points, what is a good value? ot_cost = ot.emd2(a, b, M, numItermax=2000000) - return ot_cost ** (1./order) + return (ot_cost + essential_cost) ** (1./order) diff --git a/src/python/pyproject.toml b/src/python/pyproject.toml new file mode 100644 index 00000000..a9fb4985 --- /dev/null +++ b/src/python/pyproject.toml @@ -0,0 +1,3 @@ +[build-system] +requires = ["setuptools", "wheel", "numpy>=1.15.0", "cython", "pybind11"] +build-backend = "setuptools.build_meta" diff --git a/src/python/setup.py.in b/src/python/setup.py.in index 98d058fc..23746998 100644 --- a/src/python/setup.py.in +++ b/src/python/setup.py.in @@ -41,10 +41,9 @@ for module in cython_modules: libraries=libraries, library_dirs=library_dirs, include_dirs=include_dirs, - runtime_library_dirs=runtime_library_dirs, - cython_directives = {'language_level': str(sys.version_info[0])},)) + runtime_library_dirs=runtime_library_dirs,)) -ext_modules = cythonize(ext_modules) +ext_modules = cythonize(ext_modules, compiler_directives={'language_level': str(sys.version_info[0])}) for module in pybind11_modules: my_include_dirs = include_dirs + [pybind11.get_include(False), pybind11.get_include(True)] @@ -72,7 +71,7 @@ setup( name = 'gudhi', packages=find_packages(), # find_namespace_packages(include=["gudhi*"]) author='GUDHI Editorial Board', - author_email='gudhi-contact@lists.gforge.inria.fr', + author_email='gudhi-contact@inria.fr', version='@GUDHI_VERSION@', url='https://gudhi.inria.fr/', project_urls={ @@ -83,10 +82,10 @@ setup( }, description='The Gudhi library is an open source library for ' \ 'Computational Topology and Topological Data Analysis (TDA).', + data_files=[('.', ['./introduction.rst'])], long_description_content_type='text/x-rst', long_description=long_description, ext_modules = ext_modules, - install_requires = ['numpy >= 1.9',], - setup_requires = ['cython','numpy >= 1.9','pybind11',], + install_requires = ['numpy >= 1.15.0',], package_data={"": ["*.dll"], }, ) diff --git a/src/python/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py index d0e4e9e8..29d559b3 100755 --- a/src/python/test/test_cubical_complex.py +++ b/src/python/test/test_cubical_complex.py @@ -174,3 +174,28 @@ def test_periodic_cofaces_of_persistence_pairs_when_pd_has_no_paired_birth_and_d 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])) + +def test_cubical_persistence_intervals_in_dimension(): + cub = CubicalComplex( + dimensions=[3, 3], + top_dimensional_cells=[1, 2, 3, 4, 5, 6, 7, 8, 9], + ) + cub.compute_persistence() + H0 = cub.persistence_intervals_in_dimension(0) + assert np.array_equal(H0, np.array([[ 1., float("inf")]])) + assert cub.persistence_intervals_in_dimension(1).shape == (0, 2) + +def test_periodic_cubical_persistence_intervals_in_dimension(): + cub = PeriodicCubicalComplex( + dimensions=[3, 3], + top_dimensional_cells=[1, 2, 3, 4, 5, 6, 7, 8, 9], + periodic_dimensions = [True, True] + ) + cub.compute_persistence() + H0 = cub.persistence_intervals_in_dimension(0) + assert np.array_equal(H0, np.array([[ 1., float("inf")]])) + H1 = cub.persistence_intervals_in_dimension(1) + assert np.array_equal(H1, np.array([[ 3., float("inf")], [ 7., float("inf")]])) + H2 = cub.persistence_intervals_in_dimension(2) + assert np.array_equal(H2, np.array([[ 9., float("inf")]])) + assert cub.persistence_intervals_in_dimension(3).shape == (0, 2) diff --git a/src/python/test/test_datasets_generators.py b/src/python/test/test_datasets_generators.py new file mode 100755 index 00000000..91ec4a65 --- /dev/null +++ b/src/python/test/test_datasets_generators.py @@ -0,0 +1,39 @@ +""" 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 +""" + +from gudhi.datasets.generators import points + +import pytest + +def test_sphere(): + assert points.sphere(n_samples = 10, ambient_dim = 2, radius = 1., sample = 'random').shape == (10, 2) + + with pytest.raises(ValueError): + points.sphere(n_samples = 10, ambient_dim = 2, radius = 1., sample = 'other') + +def _basic_torus(impl): + assert impl(n_samples = 64, dim = 3, sample = 'random').shape == (64, 6) + assert impl(n_samples = 64, dim = 3, sample = 'grid').shape == (64, 6) + + assert impl(n_samples = 10, dim = 4, sample = 'random').shape == (10, 8) + + # Here 1**dim < n_samples < 2**dim, the output shape is therefore (1, 2*dim) = (1, 8), where shape[0] is rounded down to the closest perfect 'dim'th power + assert impl(n_samples = 10, dim = 4, sample = 'grid').shape == (1, 8) + + with pytest.raises(ValueError): + impl(n_samples = 10, dim = 4, sample = 'other') + +def test_torus(): + for torus_impl in [points.torus, points.ctorus]: + _basic_torus(torus_impl) + # Check that the two versions (torus and ctorus) generate the same output + assert points.ctorus(n_samples = 64, dim = 3, sample = 'random').all() == points.torus(n_samples = 64, dim = 3, sample = 'random').all() + assert points.ctorus(n_samples = 64, dim = 3, sample = 'grid').all() == points.torus(n_samples = 64, dim = 3, sample = 'grid').all() + assert points.ctorus(n_samples = 10, dim = 3, sample = 'grid').all() == points.torus(n_samples = 10, dim = 3, sample = 'grid').all() diff --git a/src/python/test/test_dtm.py b/src/python/test/test_dtm.py index 0a52279e..e46d616c 100755 --- a/src/python/test/test_dtm.py +++ b/src/python/test/test_dtm.py @@ -13,6 +13,7 @@ import numpy import pytest import torch import math +import warnings def test_dtm_compare_euclidean(): @@ -87,3 +88,14 @@ def test_density(): assert density == pytest.approx(expected) density = DTMDensity(weights=[0.5, 0.5], metric="neighbors", dim=1).fit_transform(distances) assert density == pytest.approx(expected) + +def test_dtm_overflow_warnings(): + pts = numpy.array([[10., 100000000000000000000000000000.], [1000., 100000000000000000000000000.]]) + + with warnings.catch_warnings(record=True) as w: + # TODO Test "keops" implementation as well when next version of pykeops (current is 1.5) is released (should fix the problem (cf. issue #543)) + dtm = DistanceToMeasure(2, implementation="hnsw") + r = dtm.fit_transform(pts) + assert len(w) == 1 + assert issubclass(w[0].category, RuntimeWarning) + assert "Overflow" in str(w[0].message) diff --git a/src/python/test/test_reader_utils.py b/src/python/test/test_reader_utils.py index 90da6651..e96e0569 100755 --- a/src/python/test/test_reader_utils.py +++ b/src/python/test/test_reader_utils.py @@ -30,7 +30,7 @@ def test_full_square_distance_matrix_csv_file(): test_file.write("0;1;2;3;\n1;0;4;5;\n2;4;0;6;\n3;5;6;0;") test_file.close() matrix = gudhi.read_lower_triangular_matrix_from_csv_file( - csv_file="full_square_distance_matrix.csv" + csv_file="full_square_distance_matrix.csv", separator=";" ) assert matrix == [[], [1.0], [2.0, 4.0], [3.0, 5.0, 6.0]] diff --git a/src/python/test/test_representations.py b/src/python/test/test_representations.py index cda1a15b..93461f1e 100755 --- a/src/python/test/test_representations.py +++ b/src/python/test/test_representations.py @@ -3,9 +3,23 @@ import sys import matplotlib.pyplot as plt import numpy as np import pytest +import random from sklearn.cluster import KMeans +# Vectorization +from gudhi.representations import (Landscape, Silhouette, BettiCurve, ComplexPolynomial,\ + TopologicalVector, PersistenceImage, Entropy) + +# Preprocessing +from gudhi.representations import (BirthPersistenceTransform, Clamping, DiagramScaler, Padding, ProminentPoints, \ + DiagramSelector) + +# Kernel +from gudhi.representations import (PersistenceWeightedGaussianKernel, \ + PersistenceScaleSpaceKernel, SlicedWassersteinDistance,\ + SlicedWassersteinKernel, PersistenceFisherKernel, WassersteinDistance) + def test_representations_examples(): # Disable graphics for testing purposes @@ -98,3 +112,60 @@ def test_infinity(): assert c[1] == 0 assert c[7] == 3 assert c[9] == 2 + + +def test_preprocessing_empty_diagrams(): + empty_diag = np.empty(shape = [0, 2]) + assert not np.any(BirthPersistenceTransform()(empty_diag)) + assert not np.any(Clamping().fit_transform(empty_diag)) + assert not np.any(DiagramScaler()(empty_diag)) + assert not np.any(Padding()(empty_diag)) + assert not np.any(ProminentPoints()(empty_diag)) + assert not np.any(DiagramSelector()(empty_diag)) + +def pow(n): + return lambda x: np.power(x[1]-x[0],n) + +def test_vectorization_empty_diagrams(): + empty_diag = np.empty(shape = [0, 2]) + random_resolution = random.randint(50,100)*10 # between 500 and 1000 + print("resolution = ", random_resolution) + lsc = Landscape(resolution=random_resolution)(empty_diag) + assert not np.any(lsc) + assert lsc.shape[0]%random_resolution == 0 + slt = Silhouette(resolution=random_resolution, weight=pow(2))(empty_diag) + assert not np.any(slt) + assert slt.shape[0] == random_resolution + btc = BettiCurve(resolution=random_resolution)(empty_diag) + assert not np.any(btc) + assert btc.shape[0] == random_resolution + cpp = ComplexPolynomial(threshold=random_resolution, polynomial_type="T")(empty_diag) + assert not np.any(cpp) + assert cpp.shape[0] == random_resolution + tpv = TopologicalVector(threshold=random_resolution)(empty_diag) + assert tpv.shape[0] == random_resolution + assert not np.any(tpv) + prmg = PersistenceImage(resolution=[random_resolution,random_resolution])(empty_diag) + assert not np.any(prmg) + assert prmg.shape[0] == random_resolution * random_resolution + sce = Entropy(mode="scalar", resolution=random_resolution)(empty_diag) + assert not np.any(sce) + assert sce.shape[0] == 1 + scv = Entropy(mode="vector", normalized=False, resolution=random_resolution)(empty_diag) + assert not np.any(scv) + assert scv.shape[0] == random_resolution + +def test_kernel_empty_diagrams(): + empty_diag = np.empty(shape = [0, 2]) + assert SlicedWassersteinDistance(num_directions=100)(empty_diag, empty_diag) == 0. + assert SlicedWassersteinKernel(num_directions=100, bandwidth=1.)(empty_diag, empty_diag) == 1. + assert WassersteinDistance(mode="hera", delta=0.0001)(empty_diag, empty_diag) == 0. + assert WassersteinDistance(mode="pot")(empty_diag, empty_diag) == 0. + assert BottleneckDistance(epsilon=.001)(empty_diag, empty_diag) == 0. + assert BottleneckDistance()(empty_diag, empty_diag) == 0. +# PersistenceWeightedGaussianKernel(bandwidth=1., kernel_approx=None, weight=arctan(1.,1.))(empty_diag, empty_diag) +# PersistenceWeightedGaussianKernel(kernel_approx=RBFSampler(gamma=1./2, n_components=100000).fit(np.ones([1,2])), weight=arctan(1.,1.))(empty_diag, empty_diag) +# PersistenceScaleSpaceKernel(bandwidth=1.)(empty_diag, empty_diag) +# PersistenceScaleSpaceKernel(kernel_approx=RBFSampler(gamma=1./2, n_components=100000).fit(np.ones([1,2])))(empty_diag, empty_diag) +# PersistenceFisherKernel(bandwidth_fisher=1., bandwidth=1.)(empty_diag, empty_diag) +# PersistenceFisherKernel(bandwidth_fisher=1., bandwidth=1., kernel_approx=RBFSampler(gamma=1./2, n_components=100000).fit(np.ones([1,2])))(empty_diag, empty_diag) diff --git a/src/python/test/test_rips_complex.py b/src/python/test/test_rips_complex.py index b86e7498..a2f43a1b 100755 --- a/src/python/test/test_rips_complex.py +++ b/src/python/test/test_rips_complex.py @@ -133,3 +133,24 @@ def test_filtered_rips_from_distance_matrix(): assert simplex_tree.num_simplices() == 8 assert simplex_tree.num_vertices() == 4 + + +def test_sparse_with_multiplicity(): + points = [ + [3, 4], + [0.1, 2], + [0.1, 2], + [0.1, 2], + [0.1, 2], + [0.1, 2], + [0.1, 2], + [0.1, 2], + [0.1, 2], + [0.1, 2], + [0.1, 2], + [3, 4.1], + ] + rips = RipsComplex(points=points, sparse=0.01) + simplex_tree = rips.create_simplex_tree(max_dimension=2) + assert simplex_tree.num_simplices() == 7 + diag = simplex_tree.persistence() diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 92b909ca..9766ecfe 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -9,6 +9,7 @@ """ from gudhi import SimplexTree, __GUDHI_USE_EIGEN3 +import numpy as np import pytest __author__ = "Vincent Rouvreau" @@ -405,6 +406,48 @@ def test_boundaries_iterator(): with pytest.raises(RuntimeError): list(st.get_boundaries([6])) # (6) does not exist +def test_persistence_intervals_in_dimension(): + # Here is our triangulation of a 2-torus - taken from https://dioscuri-tda.org/Paris_TDA_Tutorial_2021.html + # 0-----3-----4-----0 + # | \ | \ | \ | \ | + # | \ | \ | \| \ | + # 1-----8-----7-----1 + # | \ | \ | \ | \ | + # | \ | \ | \ | \ | + # 2-----5-----6-----2 + # | \ | \ | \ | \ | + # | \ | \ | \ | \ | + # 0-----3-----4-----0 + st = SimplexTree() + st.insert([0,1,8]) + st.insert([0,3,8]) + st.insert([3,7,8]) + st.insert([3,4,7]) + st.insert([1,4,7]) + st.insert([0,1,4]) + st.insert([1,2,5]) + st.insert([1,5,8]) + st.insert([5,6,8]) + st.insert([6,7,8]) + st.insert([2,6,7]) + st.insert([1,2,7]) + st.insert([0,2,3]) + st.insert([2,3,5]) + st.insert([3,4,5]) + st.insert([4,5,6]) + st.insert([0,4,6]) + st.insert([0,2,6]) + st.compute_persistence(persistence_dim_max=True) + + H0 = st.persistence_intervals_in_dimension(0) + assert np.array_equal(H0, np.array([[ 0., float("inf")]])) + H1 = st.persistence_intervals_in_dimension(1) + assert np.array_equal(H1, np.array([[ 0., float("inf")], [ 0., float("inf")]])) + H2 = st.persistence_intervals_in_dimension(2) + assert np.array_equal(H2, np.array([[ 0., float("inf")]])) + # Test empty case + assert st.persistence_intervals_in_dimension(3).shape == (0, 2) + def test_equality_operator(): st1 = SimplexTree() st2 = SimplexTree() @@ -415,4 +458,4 @@ def test_equality_operator(): assert st1 != st2 st2.insert([1,2,3], 4.) - assert st1 == st2
\ No newline at end of file + assert st1 == st2 diff --git a/src/python/test/test_tomato.py b/src/python/test/test_tomato.py index ecab03c4..c571f799 100755 --- a/src/python/test/test_tomato.py +++ b/src/python/test/test_tomato.py @@ -37,7 +37,7 @@ def test_tomato_1(): t = Tomato(metric="euclidean", graph_type="radius", r=4.7, k=4) t.fit(a) assert t.max_weight_per_cc_.size == 2 - assert np.array_equal(t.neighbors_, [[0, 1, 2], [0, 1, 2], [0, 1, 2], [3, 4, 5, 6], [3, 4, 5], [3, 4, 5], [3, 6]]) + assert t.neighbors_ == [[0, 1, 2], [0, 1, 2], [0, 1, 2], [3, 4, 5, 6], [3, 4, 5], [3, 4, 5], [3, 6]] t.plot_diagram() t = Tomato(graph_type="radius", r=4.7, k=4, symmetrize_graph=True) diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index e3b521d6..3a004d77 100755 --- a/src/python/test/test_wasserstein_distance.py +++ b/src/python/test/test_wasserstein_distance.py @@ -5,25 +5,97 @@ Copyright (C) 2019 Inria Modification(s): + - 2020/07 Théo Lacombe: Added tests about handling essential parts in diagrams. - YYYY/MM Author: Description of the modification """ -from gudhi.wasserstein.wasserstein import _proj_on_diag +from gudhi.wasserstein.wasserstein import _proj_on_diag, _finite_part, _handle_essential_parts, _get_essential_parts +from gudhi.wasserstein.wasserstein import _warn_infty from gudhi.wasserstein import wasserstein_distance as pot from gudhi.hera import wasserstein_distance as hera import numpy as np import pytest + __author__ = "Theo Lacombe" __copyright__ = "Copyright (C) 2019 Inria" __license__ = "MIT" + def test_proj_on_diag(): dgm = np.array([[1., 1.], [1., 2.], [3., 5.]]) assert np.array_equal(_proj_on_diag(dgm), [[1., 1.], [1.5, 1.5], [4., 4.]]) empty = np.empty((0, 2)) assert np.array_equal(_proj_on_diag(empty), empty) + +def test_finite_part(): + diag = np.array([[0, 1], [3, 5], [2, np.inf], [3, np.inf], [-np.inf, 8], [-np.inf, 12], [-np.inf, -np.inf], + [np.inf, np.inf], [-np.inf, np.inf], [-np.inf, np.inf]]) + assert np.array_equal(_finite_part(diag), [[0, 1], [3, 5]]) + + +def test_handle_essential_parts(): + diag1 = np.array([[0, 1], [3, 5], + [2, np.inf], [3, np.inf], + [-np.inf, 8], [-np.inf, 12], + [-np.inf, -np.inf], + [np.inf, np.inf], + [-np.inf, np.inf], [-np.inf, np.inf]]) + + diag2 = np.array([[0, 2], [3, 5], + [2, np.inf], [4, np.inf], + [-np.inf, 8], [-np.inf, 11], + [-np.inf, -np.inf], + [np.inf, np.inf], + [-np.inf, np.inf], [-np.inf, np.inf]]) + + diag3 = np.array([[0, 2], [3, 5], + [2, np.inf], [4, np.inf], [6, np.inf], + [-np.inf, 8], [-np.inf, 11], + [-np.inf, -np.inf], + [np.inf, np.inf], + [-np.inf, np.inf], [-np.inf, np.inf]]) + + c, m = _handle_essential_parts(diag1, diag2, order=1) + assert c == pytest.approx(2, 0.0001) # Note: here c is only the cost due to essential part (thus 2, not 3) + # Similarly, the matching only corresponds to essential parts. + # Note that (-inf,-inf) and (+inf,+inf) coordinates are matched to the diagonal. + assert np.array_equal(m, [[4, 4], [5, 5], [2, 2], [3, 3], [8, 8], [9, 9], [6, -1], [7, -1], [-1, 6], [-1, 7]]) + + c, m = _handle_essential_parts(diag1, diag3, order=1) + assert c == np.inf + assert (m is None) + + +def test_get_essential_parts(): + diag1 = np.array([[0, 1], [3, 5], [2, np.inf], [3, np.inf], [-np.inf, 8], [-np.inf, 12], [-np.inf, -np.inf], + [np.inf, np.inf], [-np.inf, np.inf], [-np.inf, np.inf]]) + + diag2 = np.array([[0, 1], [3, 5], [2, np.inf], [3, np.inf]]) + + res = _get_essential_parts(diag1) + res2 = _get_essential_parts(diag2) + assert np.array_equal(res[0], [4, 5]) + assert np.array_equal(res[1], [2, 3]) + assert np.array_equal(res[2], [8, 9]) + assert np.array_equal(res[3], [6] ) + assert np.array_equal(res[4], [7] ) + + assert np.array_equal(res2[0], [] ) + assert np.array_equal(res2[1], [2, 3]) + assert np.array_equal(res2[2], [] ) + assert np.array_equal(res2[3], [] ) + assert np.array_equal(res2[4], [] ) + + +def test_warn_infty(): + assert _warn_infty(matching=False)==np.inf + c, m = _warn_infty(matching=True) + assert (c == np.inf) + assert (m is None) + + def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_matching=True): diag1 = np.array([[2.7, 3.7], [9.6, 14.0], [34.2, 34.974]]) diag2 = np.array([[2.8, 4.45], [9.5, 14.1]]) @@ -64,7 +136,7 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_mat assert wasserstein_distance(diag4, diag5) == np.inf assert wasserstein_distance(diag5, diag6, order=1, internal_p=np.inf) == approx(4.) - + assert wasserstein_distance(diag5, emptydiag) == np.inf if test_matching: match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=1., order=2)[1] @@ -78,6 +150,31 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_mat match = wasserstein_distance(diag1, diag2, matching=True, internal_p=2., order=2.)[1] assert np.array_equal(match, [[0, 0], [1, 1], [2, -1]]) + if test_matching and test_infinity: + diag7 = np.array([[0, 3], [4, np.inf], [5, np.inf]]) + diag8 = np.array([[0,1], [0, np.inf], [-np.inf, -np.inf], [np.inf, np.inf]]) + diag9 = np.array([[-np.inf, -np.inf], [np.inf, np.inf]]) + diag10 = np.array([[0,1], [-np.inf, -np.inf], [np.inf, np.inf]]) + + match = wasserstein_distance(diag5, diag6, matching=True, internal_p=2., order=2.)[1] + assert np.array_equal(match, [[0, -1], [-1,0], [-1, 1], [1, 2]]) + match = wasserstein_distance(diag5, diag7, matching=True, internal_p=2., order=2.)[1] + assert (match is None) + cost, match = wasserstein_distance(diag7, emptydiag, matching=True, internal_p=2., order=2.3) + assert (cost == np.inf) + assert (match is None) + cost, match = wasserstein_distance(emptydiag, diag7, matching=True, internal_p=2.42, order=2.) + assert (cost == np.inf) + assert (match is None) + cost, match = wasserstein_distance(diag8, diag9, matching=True, internal_p=2., order=2.) + assert (cost == np.inf) + assert (match is None) + cost, match = wasserstein_distance(diag9, diag10, matching=True, internal_p=1., order=1.) + assert (cost == 1) + assert (match == [[0, -1],[1, -1],[-1, 0], [-1, 1], [-1, 2]]) # type 4 and 5 are match to the diag anyway. + cost, match = wasserstein_distance(diag9, emptydiag, matching=True, internal_p=2., order=2.) + assert (cost == 0.) + assert (match == [[0, -1], [1, -1]]) def hera_wrap(**extra): @@ -85,15 +182,19 @@ def hera_wrap(**extra): return hera(*kargs,**kwargs,**extra) return fun + def pot_wrap(**extra): def fun(*kargs,**kwargs): return pot(*kargs,**kwargs,**extra) return fun + def test_wasserstein_distance_pot(): - _basic_wasserstein(pot, 1e-15, test_infinity=False, test_matching=True) - _basic_wasserstein(pot_wrap(enable_autodiff=True), 1e-15, test_infinity=False, test_matching=False) + _basic_wasserstein(pot, 1e-15, test_infinity=False, test_matching=True) # pot with its standard args + _basic_wasserstein(pot_wrap(enable_autodiff=True, keep_essential_parts=False), 1e-15, test_infinity=False, test_matching=False) + def test_wasserstein_distance_hera(): _basic_wasserstein(hera_wrap(delta=1e-12), 1e-12, test_matching=False) _basic_wasserstein(hera_wrap(delta=.1), .1, test_matching=False) + |
