diff options
Diffstat (limited to 'src/python')
72 files changed, 2743 insertions, 880 deletions
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 4f26481e..fbfb0d1f 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', ") @@ -103,6 +108,19 @@ if(PYTHONINTERP_FOUND) if(EAGERPY_FOUND) add_gudhi_debug_info("EagerPy version ${EAGERPY_VERSION}") endif() + 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', ") @@ -130,6 +148,10 @@ if(PYTHONINTERP_FOUND) add_gudhi_debug_info("Eigen3 version ${EIGEN3_VERSION}") # No problem, even if no CGAL found set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_EIGEN3_ENABLED', ") + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DGUDHI_USE_EIGEN3', ") + set(GUDHI_USE_EIGEN3 "True") + else (EIGEN3_FOUND) + set(GUDHI_USE_EIGEN3 "False") endif (EIGEN3_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_reader', ") @@ -144,18 +166,26 @@ 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 () if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) - set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'subsampling', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'tangential_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_witness_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_strong_witness_complex', ") endif () + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") + endif () + add_gudhi_debug_info("Boost version ${Boost_VERSION}") 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}") @@ -163,7 +193,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) @@ -171,7 +201,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. @@ -183,13 +213,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) @@ -202,7 +232,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) @@ -223,18 +253,22 @@ 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() + if(DEBUG_TRACES) + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DDEBUG_TRACES', ") + endif(DEBUG_TRACES) + if(UNIX AND WITH_GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS) set( GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}") endif(UNIX AND WITH_GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS) @@ -255,9 +289,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 @@ -267,67 +304,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)") - + # 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 "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" - ${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 "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" - ${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 5.1.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 5.1.0) + message("++ Python documentation module will not be compiled because it requires a Eigen3 and CGAL version >= 5.1.0") + set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") + endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.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") @@ -335,36 +379,36 @@ if(PYTHONINTERP_FOUND) # Test examples - if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) # Bottleneck and Alpha add_test(NAME alpha_rips_persistence_bottleneck_distance_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/alpha_rips_persistence_bottleneck_distance.py" -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -t 0.15 -d 3) - if(MATPLOTLIB_FOUND AND NUMPY_FOUND) - # Tangential - add_test(NAME tangential_complex_plain_homology_from_off_file_example_py_test - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" - ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/tangential_complex_plain_homology_from_off_file_example.py" - --no-diagram -i 2 -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off) - - add_gudhi_py_test(test_tangential_complex) - - # Witness complex AND Subsampling - add_test(NAME euclidean_strong_witness_complex_diagram_persistence_from_off_file_example_py_test - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" - ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py" - --no-diagram -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -a 1.0 -n 20 -d 2) - - add_test(NAME euclidean_witness_complex_diagram_persistence_from_off_file_example_py_test - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" - ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py" - --no-diagram -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -a 1.0 -n 20 -d 2) - endif() + endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + # Tangential + add_test(NAME tangential_complex_plain_homology_from_off_file_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/tangential_complex_plain_homology_from_off_file_example.py" + --no-diagram -i 2 -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off) + + add_gudhi_py_test(test_tangential_complex) + + # Witness complex + add_test(NAME euclidean_strong_witness_complex_diagram_persistence_from_off_file_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/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py" + --no-diagram -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -a 1.0 -n 20 -d 2) + + add_test(NAME euclidean_witness_complex_diagram_persistence_from_off_file_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/euclidean_witness_complex_diagram_persistence_from_off_file_example.py" + --no-diagram -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -a 1.0 -n 20 -d 2) # Subsampling add_gudhi_py_test(test_subsampling) @@ -374,12 +418,10 @@ if(PYTHONINTERP_FOUND) # Bottleneck add_test(NAME bottleneck_basic_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" + 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}/) @@ -387,89 +429,90 @@ if(PYTHONINTERP_FOUND) file(COPY ${CMAKE_SOURCE_DIR}/data/points/COIL_database/lucky_cat_PCA1 DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) add_test(NAME cover_complex_nerve_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/nerve_of_a_covering.py" -f human.off -c 2 -r 10 -g 0.3) add_test(NAME cover_complex_coordinate_gic_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/coordinate_graph_induced_complex.py" -f human.off -c 0 -v) add_test(NAME cover_complex_functional_gic_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/functional_graph_induced_complex.py" -o lucky_cat.off -f lucky_cat_PCA1 -v) add_test(NAME cover_complex_voronoi_gic_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/voronoi_graph_induced_complex.py" -f human.off -n 700 -v) add_gudhi_py_test(test_cover_complex) endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) - if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) # Alpha add_test(NAME alpha_complex_from_points_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${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") - if(MATPLOTLIB_FOUND AND NUMPY_FOUND) - add_test(NAME alpha_complex_diagram_persistence_from_off_file_example_py_test + add_test(NAME alpha_complex_from_generated_points_on_sphere_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" - ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/alpha_complex_diagram_persistence_from_off_file_example.py" - --no-diagram -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -a 0.6) - endif() + 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}" + ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/alpha_complex_diagram_persistence_from_off_file_example.py" + --no-diagram -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off) add_gudhi_py_test(test_alpha_complex) - endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) # 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 add_test(NAME periodic_cubical_complex_barcode_persistence_from_perseus_file_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py" --no-barcode -f ${CMAKE_SOURCE_DIR}/data/bitmap/CubicalTwoSphere.txt) - if(NUMPY_FOUND) - add_test(NAME random_cubical_complex_persistence_example_py_test - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" - ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/random_cubical_complex_persistence_example.py" - 10 10 10) - endif() + add_test(NAME random_cubical_complex_persistence_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/random_cubical_complex_persistence_example.py" + 10 10 10) add_gudhi_py_test(test_cubical_complex) # Rips - if(MATPLOTLIB_FOUND AND NUMPY_FOUND) - add_test(NAME rips_complex_diagram_persistence_from_distance_matrix_file_example_py_test - WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" - ${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) + add_test(NAME rips_complex_diagram_persistence_from_distance_matrix_file_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/rips_complex_diagram_persistence_from_distance_matrix_file_example.py" + --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} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" - ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/example/rips_complex_diagram_persistence_from_off_file_example.py - --no-diagram -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -e 0.25 -d 3) - endif() + add_test(NAME rips_complex_diagram_persistence_from_off_file_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/rips_complex_diagram_persistence_from_off_file_example.py + --no-diagram -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -e 0.25 -d 3) add_test(NAME rips_complex_from_points_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/example/rips_complex_from_points_example.py) add_gudhi_py_test(test_rips_complex) @@ -477,7 +520,7 @@ if(PYTHONINTERP_FOUND) # Simplex tree add_test(NAME simplex_tree_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/example/simplex_tree_example.py) add_gudhi_py_test(test_simplex_tree) @@ -486,7 +529,7 @@ if(PYTHONINTERP_FOUND) # Witness add_test(NAME witness_complex_from_nearest_landmark_table_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} - COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}" + COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/example/witness_complex_from_nearest_landmark_table.py) add_gudhi_py_test(test_witness_complex) @@ -495,18 +538,29 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_reader_utils) # Wasserstein - if(OT_FOUND AND PYBIND11_FOUND) - if(TORCH_FOUND AND EAGERPY_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) + + if(TORCH_FOUND AND TENSORFLOW_FOUND AND EAGERPY_FOUND) + add_gudhi_py_test(test_wasserstein_with_tensors) + endif() endif() # Representations - if(SKLEARN_FOUND AND MATPLOTLIB_FOUND) + if(SKLEARN_FOUND AND MATPLOTLIB_FOUND AND OT_FOUND AND NOT CGAL_VERSION VERSION_LESS 4.11.0) add_gudhi_py_test(test_representations) endif() + # Betti curves + if(SKLEARN_FOUND AND SCIPY_FOUND) + add_gudhi_py_test(test_betti_curve_representations) + endif() + # Time Delay add_gudhi_py_test(test_time_delay) @@ -517,7 +571,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() @@ -534,11 +588,11 @@ 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") -endif(PYTHONINTERP_FOUND) +endif(PYTHONINTERP_FOUND)
\ No newline at end of file diff --git a/src/python/doc/_templates/layout.html b/src/python/doc/_templates/layout.html index a672a281..e074b6c7 100644 --- a/src/python/doc/_templates/layout.html +++ b/src/python/doc/_templates/layout.html @@ -175,58 +175,60 @@ <h1 class="show-for-small-only"><a href="" class="icon-tree"> GUDHI library</a></h1> </li> <!-- Remove the class "menu-icon" to get rid of menu icon. Take out "Menu" to just have icon alone --> - <li class="toggle-topbar menu-icon"><a href="#"><span>Navigation</span></a></li> + <li class="toggle-topbar menu-icon"><a href="#"><span>Nav</span></a></li> </ul> <section class="top-bar-section"> <ul class="right"> <li class="divider"></li> - <li><a href="/contact/">Contact</a></li> + <li><a href="/contact/">Contact</a></li> </ul> <ul class="left"> - <li><a href="/"> <img src="/assets/img/home.png" alt=" GUDHI"> GUDHI </a></li> + <li><a href="/"> <img src="/assets/img/home.png" alt=" GUDHI"> GUDHI </a></li> <li class="divider"></li> <li class="has-dropdown"> - <a href="#">Project</a> + <a href="#">Project</a> <ul class="dropdown"> - <li><a href="/people/">People</a></li> - <li><a href="/keepintouch/">Keep in touch</a></li> - <li><a href="/partners/">Partners and Funding</a></li> - <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="/people/">People</a></li> + <li><a href="/keepintouch/">Keep in touch</a></li> + <li><a href="/partners/">Partners and Funding</a></li> + <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> <li class="has-dropdown"> - <a href="#">Download</a> + <a href="#">Download</a> <ul class="dropdown"> - <li><a href="/licensing/">Licensing</a></li> - <li><a href="https://github.com/GUDHI/gudhi-devel/releases/latest" target="_blank">Get the latest sources</a></li> - <li><a href="/conda/">Conda package</a></li> - <li><a href="/dockerfile/">Dockerfile</a></li> + <li><a href="/licensing/">Licensing</a></li> + <li><a href="https://github.com/GUDHI/gudhi-devel/releases/latest" target="_blank">Get the latest sources</a></li> + <li><a href="/conda/">Conda package</a></li> + <li><a href="https://pypi.org/project/gudhi/" target="_blank">Pip package</a></li> + <li><a href="/dockerfile/">Dockerfile</a></li> </ul> </li> <li class="divider"></li> <li class="has-dropdown"> - <a href="#">Documentation</a> + <a href="#">Documentation</a> <ul class="dropdown"> - <li><a href="/introduction/">Introduction</a></li> - <li><a href="https://gudhi.inria.fr/doc/latest/installation.html">C++ installation manual</a></li> - <li><a href="https://gudhi.inria.fr/doc/latest/">C++ documentation</a></li> - <li><a href="https://gudhi.inria.fr/python/latest/installation.html">Python installation manual</a></li> - <li><a href="https://gudhi.inria.fr/python/latest/">Python documentation</a></li> - <li><a href="/utils/">Utilities</a></li> - <li><a href="/tutorials/">Tutorials</a></li> + <li><a href="/introduction/">Introduction</a></li> + <li><a href="/doc/latest/installation.html">C++ installation manual</a></li> + <li><a href="/doc/latest/">C++ documentation</a></li> + <li><a href="/python/latest/installation.html">Python installation manual</a></li> + <li><a href="/python/latest/">Python documentation</a></li> + <li><a href="/utils/">Utilities</a></li> + <li><a href="/tutorials/">Tutorials</a></li> </ul> </li> <li class="divider"></li> - <li><a href="/interfaces/">Interfaces</a></li> + <li><a href="/interfaces/">Interfaces</a></li> <li class="divider"></li> </ul> </section> </nav> - </div><!-- /#navigation --> - <!-- GUDHI website header BEGIN --> + </div><!-- /#navigation --> + <!-- GUDHI website header END --> {%- block header %}{% endblock %} diff --git a/src/python/doc/alpha_complex_ref.rst b/src/python/doc/alpha_complex_ref.rst index 7da79543..eaa72551 100644 --- a/src/python/doc/alpha_complex_ref.rst +++ b/src/python/doc/alpha_complex_ref.rst @@ -9,6 +9,5 @@ Alpha complex reference manual .. autoclass:: gudhi.AlphaComplex :members: :undoc-members: - :show-inheritance: .. automethod:: gudhi.AlphaComplex.__init__ diff --git a/src/python/doc/alpha_complex_sum.inc b/src/python/doc/alpha_complex_sum.inc index aeab493f..5c76fd54 100644 --- a/src/python/doc/alpha_complex_sum.inc +++ b/src/python/doc/alpha_complex_sum.inc @@ -1,15 +1,15 @@ .. table:: :widths: 30 40 30 - +----------------------------------------------------------------+-------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------------+ - | .. figure:: | Alpha complex is a simplicial complex constructed from the finite | :Author: Vincent Rouvreau | - | ../../doc/Alpha_complex/alpha_complex_representation.png | cells of a Delaunay Triangulation. It has the same persistent homology | | - | :alt: Alpha complex representation | as the Čech complex and is significantly smaller. | :Since: GUDHI 2.0.0 | - | :figclass: align-center | | | - | | | :License: MIT (`GPL v3 </licensing/>`_) | - | | | | - | | | :Requires: `Eigen <installation.html#eigen>`_ :math:`\geq` 3.1.0 and `CGAL <installation.html#cgal>`_ :math:`\geq` 4.11.0 | - | | | | - +----------------------------------------------------------------+-------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------------+ - | * :doc:`alpha_complex_user` | * :doc:`alpha_complex_ref` | - +----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ + +----------------------------------------------------------------+-------------------------------------------------------------------------+------------------------------------------------------------------------------------------------------------------------+ + | .. figure:: | Alpha complex is a simplicial complex constructed from the finite | :Author: Vincent Rouvreau | + | ../../doc/Alpha_complex/alpha_complex_representation.png | cells of a Delaunay Triangulation. It has the same persistent homology | | + | :alt: Alpha complex representation | as the Čech complex and is significantly smaller. | :Since: GUDHI 2.0.0 | + | :figclass: align-center | | | + | | | :License: MIT (`GPL v3 </licensing/>`_) | + | | | | + | | | :Requires: `Eigen <installation.html#eigen>`_ :math:`\geq` 3.1.0 and `CGAL <installation.html#cgal>`_ :math:`\geq` 5.1 | + | | | | + +----------------------------------------------------------------+-------------------------------------------------------------------------+------------------------------------------------------------------------------------------------------------------------+ + | * :doc:`alpha_complex_user` | * :doc:`alpha_complex_ref` | + +----------------------------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst index fffcb3db..cfd22742 100644 --- a/src/python/doc/alpha_complex_user.rst +++ b/src/python/doc/alpha_complex_user.rst @@ -9,7 +9,7 @@ Definition .. include:: alpha_complex_sum.inc -:doc:`AlphaComplex <alpha_complex_ref>` is constructing a :doc:`SimplexTree <simplex_tree_ref>` using +:class:`~gudhi.AlphaComplex` is constructing a :doc:`SimplexTree <simplex_tree_ref>` using `Delaunay Triangulation <http://doc.cgal.org/latest/Triangulation/index.html#Chapter_Triangulations>`_ :cite:`cgal:hdj-t-19b` from the `Computational Geometry Algorithms Library <http://www.cgal.org/>`_ :cite:`cgal:eb-19b`. @@ -33,9 +33,6 @@ Remarks Using :code:`precision = 'fast'` makes the computations slightly faster, and the combinatorics are still exact, but the computation of filtration values can exceptionally be arbitrarily bad. In all cases, we still guarantee that the output is a valid filtration (faces have a filtration value no larger than their cofaces). -* For performances reasons, it is advised to use Alpha_complex with `CGAL <installation.html#cgal>`_ :math:`\geq` 5.0.0. -* The vertices in the output simplex tree are not guaranteed to match the order of the input points. One can use - :func:`~gudhi.AlphaComplex.get_point` to get the initial point back. Example from points ------------------- @@ -44,23 +41,22 @@ This example builds the alpha-complex from the given points: .. testcode:: - import gudhi - alpha_complex = gudhi.AlphaComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]]) + from gudhi import AlphaComplex + ac = AlphaComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]]) + + stree = ac.create_simplex_tree() + print('Alpha complex is of dimension ', stree.dimension(), ' - ', + stree.num_simplices(), ' simplices - ', stree.num_vertices(), ' vertices.') - 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) fmt = '%s -> %.2f' - for filtered_value in simplex_tree.get_filtration(): + for filtered_value in stree.get_filtration(): print(fmt % tuple(filtered_value)) The output is: .. testoutput:: - Alpha complex is of dimension 2 - 25 simplices - 7 vertices. + Alpha complex is of dimension 2 - 25 simplices - 7 vertices. [0] -> 0.00 [1] -> 0.00 [2] -> 0.00 @@ -163,7 +159,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 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -174,11 +173,75 @@ of speed-up, since we always first build the full filtered complex, so it is rec :paramref:`~gudhi.AlphaComplex.create_simplex_tree.max_alpha_square`. In the following example, a threshold of :math:`\alpha^2 = 32.0` is used. +Weighted version +^^^^^^^^^^^^^^^^ + +A weighted version for Alpha complex is available. It is like a usual Alpha complex, but based on a +`CGAL regular triangulation <https://doc.cgal.org/latest/Triangulation/index.html#title20>`_. + +This example builds the weighted alpha-complex of a small molecule, where atoms have different sizes. +It is taken from +`CGAL 3d weighted alpha shapes <https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title13>`_. + +Then, it is asked to display information about the alpha complex. + +.. testcode:: + + from gudhi import AlphaComplex + wgt_ac = AlphaComplex(points=[[ 1., -1., -1.], + [-1., 1., -1.], + [-1., -1., 1.], + [ 1., 1., 1.], + [ 2., 2., 2.]], + weights = [4., 4., 4., 4., 1.]) + + stree = wgt_ac.create_simplex_tree() + print('Weighted alpha complex is of dimension ', stree.dimension(), ' - ', + stree.num_simplices(), ' simplices - ', stree.num_vertices(), ' vertices.') + fmt = '%s -> %.2f' + for simplex in stree.get_simplices(): + print(fmt % tuple(simplex)) + +The output is: + +.. testoutput:: + + Weighted alpha complex is of dimension 3 - 29 simplices - 5 vertices. + [0, 1, 2, 3] -> -1.00 + [0, 1, 2] -> -1.33 + [0, 1, 3, 4] -> 95.00 + [0, 1, 3] -> -1.33 + [0, 1, 4] -> 95.00 + [0, 1] -> -2.00 + [0, 2, 3, 4] -> 95.00 + [0, 2, 3] -> -1.33 + [0, 2, 4] -> 95.00 + [0, 2] -> -2.00 + [0, 3, 4] -> 23.00 + [0, 3] -> -2.00 + [0, 4] -> 23.00 + [0] -> -4.00 + [1, 2, 3, 4] -> 95.00 + [1, 2, 3] -> -1.33 + [1, 2, 4] -> 95.00 + [1, 2] -> -2.00 + [1, 3, 4] -> 23.00 + [1, 3] -> -2.00 + [1, 4] -> 23.00 + [1] -> -4.00 + [2, 3, 4] -> 23.00 + [2, 3] -> -2.00 + [2, 4] -> 23.00 + [2] -> -4.00 + [3, 4] -> -1.00 + [3] -> -4.00 + [4] -> -1.00 Example from OFF file ^^^^^^^^^^^^^^^^^^^^^ -This example builds the alpha complex from 300 random points on a 2-torus. +This example builds the alpha complex from 300 random points on a 2-torus, given by an +`OFF file <fileformats.html#off-file-format>`_. Then, it computes the persistence diagram and displays it: @@ -186,14 +249,10 @@ Then, it computes the persistence diagram and displays it: :include-source: import matplotlib.pyplot as plt - import gudhi - alpha_complex = gudhi.AlphaComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/tore3D_300.off') - 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) - diag = simplex_tree.persistence() - gudhi.plot_persistence_diagram(diag) + import gudhi as gd + off_file = gd.__root_source_dir__ + '/data/points/tore3D_300.off' + points = gd.read_points_from_off_file(off_file = off_file) + stree = gd.AlphaComplex(points = points).create_simplex_tree() + dgm = stree.persistence() + gd.plot_persistence_diagram(dgm, legend = True) plt.show() diff --git a/src/python/doc/bottleneck_distance_user.rst b/src/python/doc/bottleneck_distance_user.rst index 6c6e08d9..7baa76cc 100644 --- a/src/python/doc/bottleneck_distance_user.rst +++ b/src/python/doc/bottleneck_distance_user.rst @@ -47,7 +47,7 @@ The following example explains how the distance is computed: :figclass: align-center The point (0, 13) is at distance 6.5 from the diagonal and more - specifically from the point (6.5, 6.5) + specifically from the point (6.5, 6.5). Basic example @@ -72,6 +72,6 @@ The output is: .. testoutput:: - Bottleneck distance approximation = 0.81 + Bottleneck distance approximation = 0.72 Bottleneck distance value = 0.75 diff --git a/src/python/doc/conf.py b/src/python/doc/conf.py index 3cc5d1d6..e69e2751 100755 --- a/src/python/doc/conf.py +++ b/src/python/doc/conf.py @@ -44,6 +44,8 @@ extensions = [ 'sphinx_paramlinks', ] +bibtex_bibfiles = ['../../biblio/bibliography.bib'] + todo_include_todos = True # plot option : do not show hyperlinks (Source code, png, hires.png, pdf) plot_html_show_source_link = False @@ -118,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/cubical_complex_user.rst b/src/python/doc/cubical_complex_user.rst index 3fd4e27a..6a211347 100644 --- a/src/python/doc/cubical_complex_user.rst +++ b/src/python/doc/cubical_complex_user.rst @@ -47,8 +47,8 @@ be a set of two elements). For further details and theory of cubical complexes, please consult :cite:`kaczynski2004computational` as well as the following paper :cite:`peikert2012topological`. -Data structure. ---------------- +Data structure +-------------- The implementation of Cubical complex provides a representation of complexes that occupy a rectangular region in :math:`\mathbb{R}^n`. This extra assumption allows for a memory efficient way of storing cubical complexes in a form @@ -77,8 +77,8 @@ Knowing the sizes of the bitmap, by a series of modulo operation, we can determi present in the product that gives the cube :math:`C`. In a similar way, we can compute boundary and the coboundary of each cube. Further details can be found in the literature. -Input Format. -------------- +Input Format +------------ In the current implantation, filtration is given at the maximal cubes, and it is then extended by the lower star filtration to all cubes. There are a number of constructors that can be used to construct cubical complex by users @@ -108,8 +108,8 @@ the program output is: Cubical complex is of dimension 2 - 49 simplices. -Periodic boundary conditions. ------------------------------ +Periodic boundary conditions +---------------------------- Often one would like to impose periodic boundary conditions to the cubical complex (cf. :doc:`periodic_cubical_complex_ref`). @@ -154,7 +154,13 @@ the program output is: Periodic cubical complex is of dimension 2 - 42 simplices. -Examples. ---------- +Examples +-------- End user programs are available in python/example/ folder. + +Tutorial +-------- + +This `notebook <https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-cubical-complexes.ipynb>`_ +explains how to represent sublevels sets of functions using cubical complexes.
\ No newline at end of file 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 a42227e3..1442f185 100644 --- a/src/python/doc/examples.rst +++ b/src/python/doc/examples.rst @@ -7,27 +7,31 @@ Examples .. only:: builder_html - * :download:`rips_complex_from_points_example.py <../example/rips_complex_from_points_example.py>` + * :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:`simplex_tree_example.py <../example/simplex_tree_example.py>` * :download:`alpha_rips_persistence_bottleneck_distance.py <../example/alpha_rips_persistence_bottleneck_distance.py>` - * :download:`tangential_complex_plain_homology_from_off_file_example.py <../example/tangential_complex_plain_homology_from_off_file_example.py>` - * :download:`alpha_complex_diagram_persistence_from_off_file_example.py <../example/alpha_complex_diagram_persistence_from_off_file_example.py>` - * :download:`periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py <../example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py>` * :download:`bottleneck_basic_example.py <../example/bottleneck_basic_example.py>` - * :download:`gudhi_graphical_tools_example.py <../example/gudhi_graphical_tools_example.py>` - * :download:`plot_simplex_tree_dim012.py <../example/plot_simplex_tree_dim012.py>` - * :download:`plot_rips_complex.py <../example/plot_rips_complex.py>` - * :download:`plot_alpha_complex.py <../example/plot_alpha_complex.py>` - * :download:`witness_complex_from_nearest_landmark_table.py <../example/witness_complex_from_nearest_landmark_table.py>` + * :download:`coordinate_graph_induced_complex.py <../example/coordinate_graph_induced_complex.py>` + * :download:`diagram_vectorizations_distances_kernels.py <../example/diagram_vectorizations_distances_kernels.py>` * :download:`euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py <../example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py>` * :download:`euclidean_witness_complex_diagram_persistence_from_off_file_example.py <../example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py>` - * :download:`rips_complex_diagram_persistence_from_off_file_example.py <../example/rips_complex_diagram_persistence_from_off_file_example.py>` + * :download:`functional_graph_induced_complex.py <../example/functional_graph_induced_complex.py>` + * :download:`gudhi_graphical_tools_example.py <../example/gudhi_graphical_tools_example.py>` + * :download:`nerve_of_a_covering.py <../example/nerve_of_a_covering.py>` + * :download:`periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py <../example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py>` + * :download:`plot_alpha_complex.py <../example/plot_alpha_complex.py>` + * :download:`plot_rips_complex.py <../example/plot_rips_complex.py>` + * :download:`plot_simplex_tree_dim012.py <../example/plot_simplex_tree_dim012.py>` + * :download:`random_cubical_complex_persistence_example.py <../example/random_cubical_complex_persistence_example.py>` + * :download:`rips_complex_diagram_persistence_from_correlation_matrix_file_example.py <../example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py>` * :download:`rips_complex_diagram_persistence_from_distance_matrix_file_example.py <../example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py>` + * :download:`rips_complex_diagram_persistence_from_off_file_example.py <../example/rips_complex_diagram_persistence_from_off_file_example.py>` + * :download:`rips_complex_edge_collapse_example.py <../example/rips_complex_edge_collapse_example.py>` + * :download:`rips_complex_from_points_example.py <../example/rips_complex_from_points_example.py>` * :download:`rips_persistence_diagram.py <../example/rips_persistence_diagram.py>` + * :download:`simplex_tree_example.py <../example/simplex_tree_example.py>` * :download:`sparse_rips_persistence_diagram.py <../example/sparse_rips_persistence_diagram.py>` - * :download:`random_cubical_complex_persistence_example.py <../example/random_cubical_complex_persistence_example.py>` - * :download:`coordinate_graph_induced_complex.py <../example/coordinate_graph_induced_complex.py>` - * :download:`functional_graph_induced_complex.py <../example/functional_graph_induced_complex.py>` + * :download:`tangential_complex_plain_homology_from_off_file_example.py <../example/tangential_complex_plain_homology_from_off_file_example.py>` * :download:`voronoi_graph_induced_complex.py <../example/voronoi_graph_induced_complex.py>` - * :download:`nerve_of_a_covering.py <../example/nerve_of_a_covering.py>` + * :download:`witness_complex_from_nearest_landmark_table.py <../example/witness_complex_from_nearest_landmark_table.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 525ca84e..cff84691 100644 --- a/src/python/doc/installation.rst +++ b/src/python/doc/installation.rst @@ -33,25 +33,19 @@ Compiling These instructions are for people who want to compile gudhi from source, they are unnecessary if you installed a binary package of Gudhi as above. They assume that you have downloaded a `release <https://github.com/GUDHI/gudhi-devel/releases>`_, -with a name like `gudhi.3.2.0.tar.gz`, then run `tar xf gudhi.3.2.0.tar.gz`, which -created a directory `gudhi.3.2.0`, hereinafter referred to as `/path-to-gudhi/`. +with a name like `gudhi.3.X.Y.tar.gz`, then run `tar xf gudhi.3.X.Y.tar.gz`, which +created a directory `gudhi.3.X.Y`, hereinafter referred to as `/path-to-gudhi/`. If you are instead using a git checkout, beware that the paths are a bit different, and in particular the `python/` subdirectory is actually `src/python/` 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.1 to generate makefiles, -`NumPy <http://numpy.org>`_, `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 -Studio 2017 or later. +The library uses c++14 and requires `Boost <https://www.boost.org/>`_ :math:`\geq` 1.66.0, +`CMake <https://www.cmake.org/>`_ :math:`\geq` 3.5 to generate makefiles, +Python :math:`\geq` 3.5, `NumPy <http://numpy.org>`_ :math:`\geq` 1.15.0, `Cython <https://www.cython.org/>`_ +:math:`\geq` 0.27 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 Studio 2017 or later. -On `Windows <https://wiki.python.org/moin/WindowsCompilers>`_ , only Python -:math:`\geq` 3.5 are available because of the required Visual Studio version. - -On other systems, if you have several Python/python installed, the version 2.X -will be used by default, but you can force it by adding +If you have several Python/python installed, the version 2.X may be used by default, but you can force it by adding :code:`-DPython_ADDITIONAL_VERSIONS=3` to the cmake command. GUDHI Python module compilation @@ -65,7 +59,7 @@ one can build the GUDHI Python module, by running the following commands in a te cd /path-to-gudhi/ mkdir build cd build/ - cmake .. + cmake -DCMAKE_BUILD_TYPE=Release .. cd python make @@ -99,20 +93,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 =========== @@ -148,60 +136,71 @@ If :code:`import gudhi` succeeds, please have a look to debug information: .. code-block:: python - import gudhi - print(gudhi.__debug_info__) + import gudhi as gd + print(gd.__debug_info__) + print("+ Installed modules are: " + gd.__available_modules) + print("+ Missing modules are: " + gd.__missing_modules) You shall have something like: .. code-block:: none - Python version 2.7.15 - Cython version 0.26.1 - Numpy version 1.14.1 - Eigen3 version 3.1.1 - Installed modules are: off_reader;simplex_tree;rips_complex; - cubical_complex;periodic_cubical_complex;reader_utils;witness_complex; - strong_witness_complex;alpha_complex; - Missing modules are: bottleneck_distance;nerve_gic;subsampling; - tangential_complex;persistence_graphical_tools; - euclidean_witness_complex;euclidean_strong_witness_complex; - CGAL version 4.7.1000 - GMP_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmp.so - GMPXX_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmpxx.so - TBB version 9107 found and used + Pybind11 version 2.8.1 + Python version 3.7.12 + Cython version 0.29.25 + Numpy version 1.21.4 + Boost version 1.77.0 + + Installed modules are: off_reader;simplex_tree;rips_complex;cubical_complex;periodic_cubical_complex; + persistence_graphical_tools;reader_utils;witness_complex;strong_witness_complex; + + Missing modules are: bottleneck;nerve_gic;subsampling;tangential_complex;alpha_complex;euclidean_witness_complex; + euclidean_strong_witness_complex; -Here, you can see that bottleneck_distance, nerve_gic, subsampling and -tangential_complex are missing because of the CGAL version. -persistence_graphical_tools is not available as matplotlib is not -available. +Here, you can see that the modules that need CGAL are missing, because CGAL is not installed. +:code:`persistence_graphical_tools` is installed, but +`its functions <https://gudhi.inria.fr/python/latest/persistence_graphical_tools_ref.html>`_ will produce an error as +matplotlib is not available. Unitary tests cannot be run as pytest is missing. A complete configuration would be : .. code-block:: none - Python version 3.6.5 - Cython version 0.28.2 - Pytest version 3.3.2 - Matplotlib version 2.2.2 - Numpy version 1.14.5 - Eigen3 version 3.3.4 - Installed modules are: off_reader;simplex_tree;rips_complex; - cubical_complex;periodic_cubical_complex;persistence_graphical_tools; - reader_utils;witness_complex;strong_witness_complex; - persistence_graphical_tools;bottleneck_distance;nerve_gic;subsampling; - tangential_complex;alpha_complex;euclidean_witness_complex; - euclidean_strong_witness_complex; - CGAL header only version 4.11.0 + Pybind11 version 2.8.1 + Python version 3.9.7 + Cython version 0.29.24 + Pytest version 6.2.5 + Matplotlib version 3.5.0 + Numpy version 1.21.4 + Scipy version 1.7.3 + Scikit-learn version 1.0.1 + POT version 0.8.0 + HNSWlib found + PyKeOps version [pyKeOps]: 1.5 + EagerPy version 0.30.0 + TensorFlow version 2.7.0 + Sphinx version 4.3.0 + Sphinx-paramlinks version 0.5.2 + python_docs_theme found + Eigen3 version 3.4.0 + Boost version 1.74.0 + CGAL version 5.3 GMP_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmp.so GMPXX_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmpxx.so + MPFR_LIBRARIES = /usr/lib/x86_64-linux-gnu/libmpfr.so TBB version 9107 found and used + + Installed modules are: bottleneck;off_reader;simplex_tree;rips_complex;cubical_complex;periodic_cubical_complex; + persistence_graphical_tools;reader_utils;witness_complex;strong_witness_complex;nerve_gic;subsampling; + tangential_complex;alpha_complex;euclidean_witness_complex;euclidean_strong_witness_complex; + + Missing modules are: + 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. @@ -323,6 +322,35 @@ The following examples require the `Matplotlib <http://matplotlib.org>`_: * :download:`euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py <../example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py>` * :download:`euclidean_witness_complex_diagram_persistence_from_off_file_example.py <../example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py>` +LaTeX +~~~~~ + +If a sufficiently complete LaTeX toolchain is available (including dvipng and ghostscript), the LaTeX option of +matplotlib is enabled for prettier captions (cf. +`matplotlib text rendering with LaTeX <https://matplotlib.org/3.3.0/tutorials/text/usetex.html>`_). +It also requires `type1cm` LaTeX package (not detected by matplotlib). + +If you are facing issues with LaTeX rendering, like this one: + +.. code-block:: none + + Traceback (most recent call last): + File "/usr/lib/python3/dist-packages/matplotlib/texmanager.py", line 302, in _run_checked_subprocess + report = subprocess.check_output(command, + ... + ! LaTeX Error: File `type1cm.sty' not found. + ... + +This is because the LaTeX package is not installed on your system. On Ubuntu systems you can install texlive-full +(for all LaTeX packages), or more specific packages like texlive-latex-extra, cm-super. + +You can still deactivate LaTeX rendering by saying: + +.. code-block:: python + + import gudhi as gd + gd.persistence_graphical_tools._gudhi_matplotlib_use_tex=False + PyKeOps ------- @@ -334,7 +362,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 @@ -365,11 +393,17 @@ mathematics, science, and engineering. :class:`~gudhi.point_cloud.knn.KNearestNeighbors` can use the Python package `SciPy <http://scipy.org>`_ as a backend if explicitly requested. +TensorFlow +---------- + +`TensorFlow <https://www.tensorflow.org>`_ is currently only used in some automatic differentiation tests. + 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/persistence_graphical_tools_user.rst b/src/python/doc/persistence_graphical_tools_user.rst index b5a38eb1..d95b9d2b 100644 --- a/src/python/doc/persistence_graphical_tools_user.rst +++ b/src/python/doc/persistence_graphical_tools_user.rst @@ -90,3 +90,14 @@ If you want more information on a specific dimension, for instance: gudhi.plot_persistence_density(persistence=pers_diag, dimension=1, legend=True, axes=axes[1]) plt.show() + +LaTeX support +------------- + +If you are facing issues with `LaTeX <installation.html#latex>`_ rendering, you can still deactivate LaTeX rendering by +saying: + +.. code-block:: python + + import gudhi + gudhi.persistence_graphical_tools._gudhi_matplotlib_use_tex=False diff --git a/src/python/doc/representations.rst b/src/python/doc/representations.rst index 041e3247..b0477197 100644 --- a/src/python/doc/representations.rst +++ b/src/python/doc/representations.rst @@ -12,11 +12,45 @@ This module, originally available at https://github.com/MathieuCarriere/sklearn- A diagram is represented as a numpy array of shape (n,2), as can be obtained from :func:`~gudhi.SimplexTree.persistence_intervals_in_dimension` for instance. Points at infinity are represented as a numpy array of shape (n,1), storing only the birth time. The classes in this module can handle several persistence diagrams at once. In that case, the diagrams are provided as a list of numpy arrays. Note that it is not necessary for the diagrams to have the same number of points, i.e., for the corresponding arrays to have the same number of rows: all classes can handle arrays with different shapes. -A small example is provided +Examples +-------- -.. only:: builder_html +Landscapes +^^^^^^^^^^ - * :download:`diagram_vectorizations_distances_kernels.py <../example/diagram_vectorizations_distances_kernels.py>` +This example computes the first two Landscapes associated to a persistence diagram with four points. The landscapes are evaluated on ten samples, leading to two vectors with ten coordinates each, that are eventually concatenated in order to produce a single vector representation. + +.. testcode:: + + import numpy as np + from gudhi.representations import Landscape + # A single diagram with 4 points + D = np.array([[0.,4.],[1.,2.],[3.,8.],[6.,8.]]) + diags = [D] + l=Landscape(num_landscapes=2,resolution=10).fit_transform(diags) + print(l) + +The output is: + +.. testoutput:: + + [[1.02851895 2.05703791 2.57129739 1.54277843 0.89995409 1.92847304 + 2.95699199 3.08555686 2.05703791 1.02851895 0. 0.64282435 + 0. 0. 0.51425948 0. 0. 0. + 0.77138922 1.02851895]] + +Various kernels +^^^^^^^^^^^^^^^ + +This small example is also provided +:download:`diagram_vectorizations_distances_kernels.py <../example/diagram_vectorizations_distances_kernels.py>` + +Machine Learning and Topological Data Analysis +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This `notebook <https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-representations.ipynb>`_ explains how to +efficiently combine machine learning and topological data analysis with the +:doc:`representations module<representations>`. Preprocessing @@ -46,27 +80,3 @@ Metrics :members: :special-members: :show-inheritance: - -Basic example -------------- - -This example computes the first two Landscapes associated to a persistence diagram with four points. The landscapes are evaluated on ten samples, leading to two vectors with ten coordinates each, that are eventually concatenated in order to produce a single vector representation. - -.. testcode:: - - import numpy as np - from gudhi.representations import Landscape - # A single diagram with 4 points - D = np.array([[0.,4.],[1.,2.],[3.,8.],[6.,8.]]) - diags = [D] - l=Landscape(num_landscapes=2,resolution=10).fit_transform(diags) - print(l) - -The output is: - -.. testoutput:: - - [[1.02851895 2.05703791 2.57129739 1.54277843 0.89995409 1.92847304 - 2.95699199 3.08555686 2.05703791 1.02851895 0. 0.64282435 - 0. 0. 0.51425948 0. 0. 0. - 0.77138922 1.02851895]] diff --git a/src/python/doc/representations_sum.inc b/src/python/doc/representations_sum.inc index 323a0920..4298aea9 100644 --- a/src/python/doc/representations_sum.inc +++ b/src/python/doc/representations_sum.inc @@ -2,7 +2,7 @@ :widths: 30 40 30 +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------+ - | .. figure:: | Vectorizations, distances and kernels that work on persistence | :Author: Mathieu Carrière | + | .. figure:: | Vectorizations, distances and kernels that work on persistence | :Author: Mathieu Carrière, Martin Royer | | img/sklearn-tda.png | diagrams, compatible with scikit-learn. | | | | | :Since: GUDHI 3.1.0 | | | | | diff --git a/src/python/doc/rips_complex_sum.inc b/src/python/doc/rips_complex_sum.inc index c123ea2a..2cb24990 100644 --- a/src/python/doc/rips_complex_sum.inc +++ b/src/python/doc/rips_complex_sum.inc @@ -1,14 +1,14 @@ .. table:: :widths: 30 40 30 - +----------------------------------------------------------------+------------------------------------------------------------------------+----------------------------------------------------------------------+ - | .. figure:: | The Vietoris-Rips complex is a simplicial complex built as the | :Authors: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse | - | ../../doc/Rips_complex/rips_complex_representation.png | clique-complex of a proximity graph. | | - | :figclass: align-center | | :Since: GUDHI 2.0.0 | - | | We also provide sparse approximations, to speed-up the computation | | - | | of persistent homology, and weighted versions, which are more robust | :License: MIT | - | | to outliers. | | - | | | | - +----------------------------------------------------------------+------------------------------------------------------------------------+----------------------------------------------------------------------+ - | * :doc:`rips_complex_user` | * :doc:`rips_complex_ref` | - +----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------+ + +----------------------------------------------------------------+------------------------------------------------------------------------+----------------------------------------------------------------------------------+ + | .. figure:: | The Vietoris-Rips complex is a simplicial complex built as the | :Authors: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse, Yuichi Ike | + | ../../doc/Rips_complex/rips_complex_representation.png | clique-complex of a proximity graph. | | + | :figclass: align-center | | :Since: GUDHI 2.0.0 | + | | We also provide sparse approximations, to speed-up the computation | | + | | of persistent homology, and weighted versions, which are more robust | :License: MIT | + | | to outliers. | | + | | | | + +----------------------------------------------------------------+------------------------------------------------------------------------+----------------------------------------------------------------------------------+ + | * :doc:`rips_complex_user` | * :doc:`rips_complex_ref` | + +----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/doc/rips_complex_user.rst b/src/python/doc/rips_complex_user.rst index 6048cc4e..27d218d4 100644 --- a/src/python/doc/rips_complex_user.rst +++ b/src/python/doc/rips_complex_user.rst @@ -7,9 +7,9 @@ Rips complex user manual Definition ---------- -==================================================================== ================================ ====================== -:Authors: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse :Since: GUDHI 2.0.0 :License: GPL v3 -==================================================================== ================================ ====================== +================================================================================ ================================ ====================== +:Authors: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse, Yuichi Ike :Since: GUDHI 2.0.0 :License: GPL v3 +================================================================================ ================================ ====================== +-------------------------------------------+----------------------------------------------------------------------+ | :doc:`rips_complex_user` | :doc:`rips_complex_ref` | diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst index 96ec7872..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 ----------- @@ -175,3 +192,10 @@ The output is: [[0.27916667 0.55416667] [0.7375 0.7625 ] [0.2375 0.2625 ]] + +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. diff --git a/src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py b/src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py index 727af4fa..c96121a6 100755 --- a/src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py +++ b/src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py @@ -1,10 +1,7 @@ #!/usr/bin/env python import argparse -import errno -import os -import matplotlib.pyplot as plot -import gudhi +import gudhi as gd """ This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. @@ -26,12 +23,12 @@ parser = argparse.ArgumentParser( description="AlphaComplex creation from " "points read in a OFF file.", epilog="Example: " "example/alpha_complex_diagram_persistence_from_off_file_example.py " - "-f ../data/points/tore3D_300.off -a 0.6" + "-f ../data/points/tore3D_300.off" "- Constructs a alpha complex with the " "points from the given OFF file.", ) parser.add_argument("-f", "--file", type=str, required=True) -parser.add_argument("-a", "--max_alpha_square", type=float, default=0.5) +parser.add_argument("-a", "--max_alpha_square", type=float, required=False) parser.add_argument("-b", "--band", type=float, default=0.0) parser.add_argument( "--no-diagram", @@ -42,33 +39,24 @@ parser.add_argument( args = parser.parse_args() -with open(args.file, "r") as f: - first_line = f.readline() - if (first_line == "OFF\n") or (first_line == "nOFF\n"): - print("##############################################################") - print("AlphaComplex creation from points read in a OFF file") - - message = "AlphaComplex with max_edge_length=" + repr(args.max_alpha_square) - print(message) - - alpha_complex = gudhi.AlphaComplex(off_file=args.file) - simplex_tree = alpha_complex.create_simplex_tree( - max_alpha_square=args.max_alpha_square - ) - - message = "Number of simplices=" + repr(simplex_tree.num_simplices()) - print(message) - - diag = simplex_tree.persistence() - - print("betti_numbers()=") - print(simplex_tree.betti_numbers()) - - if args.no_diagram == False: - gudhi.plot_persistence_diagram(diag, band=args.band) - plot.show() - else: - raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), - args.file) - - f.close() +print("##############################################################") +print("AlphaComplex creation from points read in a OFF file") + +points = gd.read_points_from_off_file(off_file = args.file) +alpha_complex = gd.AlphaComplex(points = points) +if args.max_alpha_square is not None: + print("with max_edge_length=", args.max_alpha_square) + simplex_tree = alpha_complex.create_simplex_tree( + max_alpha_square=args.max_alpha_square + ) +else: + simplex_tree = alpha_complex.create_simplex_tree() + +print("Number of simplices=", simplex_tree.num_simplices()) + +diag = simplex_tree.persistence() +print("betti_numbers()=", simplex_tree.betti_numbers()) +if args.no_diagram == False: + import matplotlib.pyplot as plot + gd.plot_persistence_diagram(diag, band=args.band) + plot.show() 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/alpha_complex_from_points_example.py b/src/python/example/alpha_complex_from_points_example.py index 465632eb..5d5ca66a 100755 --- a/src/python/example/alpha_complex_from_points_example.py +++ b/src/python/example/alpha_complex_from_points_example.py @@ -19,7 +19,7 @@ __license__ = "MIT" print("#####################################################################") print("AlphaComplex creation from points") alpha_complex = AlphaComplex(points=[[0, 0], [1, 0], [0, 1], [1, 1]]) -simplex_tree = alpha_complex.create_simplex_tree(max_alpha_square=60.0) +simplex_tree = alpha_complex.create_simplex_tree() if simplex_tree.find([0, 1]): print("[0, 1] Found !!") diff --git a/src/python/example/alpha_rips_persistence_bottleneck_distance.py b/src/python/example/alpha_rips_persistence_bottleneck_distance.py index 3e12b0d5..6b97fb3b 100755 --- a/src/python/example/alpha_rips_persistence_bottleneck_distance.py +++ b/src/python/example/alpha_rips_persistence_bottleneck_distance.py @@ -1,10 +1,8 @@ #!/usr/bin/env python -import gudhi +import gudhi as gd import argparse import math -import errno -import os import numpy as np """ This file is part of the Gudhi Library - https://gudhi.inria.fr/ - @@ -37,70 +35,60 @@ parser.add_argument("-t", "--threshold", type=float, default=0.5) parser.add_argument("-d", "--max_dimension", type=int, default=1) args = parser.parse_args() -with open(args.file, "r") as f: - first_line = f.readline() - if (first_line == "OFF\n") or (first_line == "nOFF\n"): - point_cloud = gudhi.read_points_from_off_file(off_file=args.file) - print("##############################################################") - print("RipsComplex creation from points read in a OFF file") +point_cloud = gd.read_points_from_off_file(off_file=args.file) +print("##############################################################") +print("RipsComplex creation from points read in a OFF file") - message = "RipsComplex with max_edge_length=" + repr(args.threshold) - print(message) +message = "RipsComplex with max_edge_length=" + repr(args.threshold) +print(message) - rips_complex = gudhi.RipsComplex( - points=point_cloud, max_edge_length=args.threshold - ) - - rips_stree = rips_complex.create_simplex_tree( - max_dimension=args.max_dimension) - - message = "Number of simplices=" + repr(rips_stree.num_simplices()) - print(message) - - rips_stree.compute_persistence() - - print("##############################################################") - print("AlphaComplex creation from points read in a OFF file") - - message = "AlphaComplex with max_edge_length=" + repr(args.threshold) - print(message) - - alpha_complex = gudhi.AlphaComplex(points=point_cloud) - alpha_stree = alpha_complex.create_simplex_tree( - max_alpha_square=(args.threshold * args.threshold) - ) - - message = "Number of simplices=" + repr(alpha_stree.num_simplices()) - print(message) +rips_complex = gd.RipsComplex( + points=point_cloud, max_edge_length=args.threshold +) - alpha_stree.compute_persistence() +rips_stree = rips_complex.create_simplex_tree( + max_dimension=args.max_dimension) - max_b_distance = 0.0 - for dim in range(args.max_dimension): - # Alpha persistence values needs to be transform because filtration - # values are alpha square values - alpha_intervals = np.sqrt(alpha_stree.persistence_intervals_in_dimension(dim)) +message = "Number of simplices=" + repr(rips_stree.num_simplices()) +print(message) - rips_intervals = rips_stree.persistence_intervals_in_dimension(dim) - bottleneck_distance = gudhi.bottleneck_distance( - rips_intervals, alpha_intervals - ) - message = ( - "In dimension " - + repr(dim) - + ", bottleneck distance = " - + repr(bottleneck_distance) - ) - print(message) - max_b_distance = max(bottleneck_distance, max_b_distance) +rips_stree.compute_persistence() - print("==============================================================") - message = "Bottleneck distance is " + repr(max_b_distance) - print(message) +print("##############################################################") +print("AlphaComplex creation from points read in a OFF file") - else: - raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), - args.file) +message = "AlphaComplex with max_edge_length=" + repr(args.threshold) +print(message) +alpha_complex = gd.AlphaComplex(points=point_cloud) +alpha_stree = alpha_complex.create_simplex_tree( + max_alpha_square=(args.threshold * args.threshold) +) - f.close() +message = "Number of simplices=" + repr(alpha_stree.num_simplices()) +print(message) + +alpha_stree.compute_persistence() + +max_b_distance = 0.0 +for dim in range(args.max_dimension): + # Alpha persistence values needs to be transform because filtration + # values are alpha square values + alpha_intervals = np.sqrt(alpha_stree.persistence_intervals_in_dimension(dim)) + + rips_intervals = rips_stree.persistence_intervals_in_dimension(dim) + bottleneck_distance = gd.bottleneck_distance( + rips_intervals, alpha_intervals + ) + message = ( + "In dimension " + + repr(dim) + + ", bottleneck distance = " + + repr(bottleneck_distance) + ) + print(message) + max_b_distance = max(bottleneck_distance, max_b_distance) + +print("==============================================================") +message = "Bottleneck distance is " + repr(max_b_distance) +print(message) diff --git a/src/python/example/diagram_vectorizations_distances_kernels.py b/src/python/example/diagram_vectorizations_distances_kernels.py index c4a71a7a..2801576e 100755 --- a/src/python/example/diagram_vectorizations_distances_kernels.py +++ b/src/python/example/diagram_vectorizations_distances_kernels.py @@ -5,11 +5,11 @@ import numpy as np from sklearn.kernel_approximation import RBFSampler from sklearn.preprocessing import MinMaxScaler -from gudhi.representations import DiagramSelector, Clamping, Landscape, Silhouette, BettiCurve, ComplexPolynomial,\ +from gudhi.representations import (DiagramSelector, Clamping, Landscape, Silhouette, BettiCurve, ComplexPolynomial,\ TopologicalVector, DiagramScaler, BirthPersistenceTransform,\ PersistenceImage, PersistenceWeightedGaussianKernel, Entropy, \ PersistenceScaleSpaceKernel, SlicedWassersteinDistance,\ - SlicedWassersteinKernel, BottleneckDistance, PersistenceFisherKernel, WassersteinDistance + SlicedWassersteinKernel, PersistenceFisherKernel, WassersteinDistance) D1 = np.array([[0.,4.],[1.,2.],[3.,8.],[6.,8.], [0., np.inf], [5., np.inf]]) @@ -93,14 +93,21 @@ print("SW distance is " + str(sW(D1, D2))) SW = SlicedWassersteinKernel(num_directions=100, bandwidth=1.) print("SW kernel is " + str(SW(D1, D2))) -W = WassersteinDistance(order=2, internal_p=2, mode="pot") -print("Wasserstein distance (POT) is " + str(W(D1, D2))) +try: + W = WassersteinDistance(order=2, internal_p=2, mode="pot") + print("Wasserstein distance (POT) is " + str(W(D1, D2))) +except ImportError: + print("WassersteinDistance (POT) is not available, you may be missing pot.") W = WassersteinDistance(order=2, internal_p=2, mode="hera", delta=0.0001) print("Wasserstein distance (hera) is " + str(W(D1, D2))) -W = BottleneckDistance(epsilon=.001) -print("Bottleneck distance is " + str(W(D1, D2))) +try: + from gudhi.representations import BottleneckDistance + W = BottleneckDistance(epsilon=.001) + print("Bottleneck distance is " + str(W(D1, D2))) +except ImportError: + print("BottleneckDistance is not available, you may be missing CGAL.") PF = PersistenceFisherKernel(bandwidth_fisher=1., bandwidth=1.) print("PF kernel is " + str(PF(D1, D2))) diff --git a/src/python/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py b/src/python/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py index e1e572df..4e97cfe3 100755 --- a/src/python/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py +++ b/src/python/example/euclidean_strong_witness_complex_diagram_persistence_from_off_file_example.py @@ -3,7 +3,6 @@ import argparse import errno import os -import matplotlib.pyplot as plot import gudhi """ This file is part of the Gudhi Library - https://gudhi.inria.fr/ - @@ -82,6 +81,7 @@ with open(args.file, "r") as f: print(simplex_tree.betti_numbers()) if args.no_diagram == False: + import matplotlib.pyplot as plot gudhi.plot_persistence_diagram(diag, band=args.band) plot.show() else: diff --git a/src/python/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py b/src/python/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py index 58cb2bb5..29076c74 100755 --- a/src/python/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py +++ b/src/python/example/euclidean_witness_complex_diagram_persistence_from_off_file_example.py @@ -3,7 +3,6 @@ import argparse import errno import os -import matplotlib.pyplot as plot import gudhi """ This file is part of the Gudhi Library - https://gudhi.inria.fr/ - @@ -79,6 +78,7 @@ with open(args.file, "r") as f: print(simplex_tree.betti_numbers()) if args.no_diagram == False: + import matplotlib.pyplot as plot gudhi.plot_persistence_diagram(diag, band=args.band) plot.show() else: diff --git a/src/python/example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py b/src/python/example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py index 499171df..ee3290c6 100755 --- a/src/python/example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py +++ b/src/python/example/periodic_cubical_complex_barcode_persistence_from_perseus_file_example.py @@ -1,7 +1,6 @@ #!/usr/bin/env python import argparse -import matplotlib.pyplot as plot import errno import os import gudhi @@ -75,6 +74,7 @@ if is_file_perseus(args.file): print("betti_numbers()=") print(periodic_cubical_complex.betti_numbers()) if args.no_barcode == False: + import matplotlib.pyplot as plot gudhi.plot_persistence_barcode(diag) plot.show() else: diff --git a/src/python/example/plot_alpha_complex.py b/src/python/example/plot_alpha_complex.py index 99c18a7c..0924619b 100755 --- a/src/python/example/plot_alpha_complex.py +++ b/src/python/example/plot_alpha_complex.py @@ -1,8 +1,9 @@ #!/usr/bin/env python import numpy as np -import gudhi -ac = gudhi.AlphaComplex(off_file='../../data/points/tore3D_1307.off') +import gudhi as gd +points = gd.read_points_from_off_file(off_file = '../../data/points/tore3D_1307.off') +ac = gd.AlphaComplex(points = points) st = ac.create_simplex_tree() points = np.array([ac.get_point(i) for i in range(st.num_vertices())]) # We want to plot the alpha-complex with alpha=0.1. diff --git a/src/python/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py b/src/python/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py index 1acb187c..ea2eb7e1 100755 --- a/src/python/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py +++ b/src/python/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py @@ -2,7 +2,6 @@ import sys import argparse -import matplotlib.pyplot as plot import gudhi """ This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. @@ -84,5 +83,6 @@ invert_diag = [ ] if args.no_diagram == False: + import matplotlib.pyplot as plot gudhi.plot_persistence_diagram(invert_diag, band=args.band) plot.show() 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 79ccca96..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 @@ -1,7 +1,6 @@ #!/usr/bin/env python import argparse -import matplotlib.pyplot as plot import gudhi """ This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. @@ -22,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) @@ -45,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 ) @@ -60,5 +60,6 @@ print("betti_numbers()=") print(simplex_tree.betti_numbers()) if args.no_diagram == False: + import matplotlib.pyplot as plot gudhi.plot_persistence_diagram(diag, band=args.band) plot.show() diff --git a/src/python/example/rips_complex_diagram_persistence_from_off_file_example.py b/src/python/example/rips_complex_diagram_persistence_from_off_file_example.py index 6f992508..e80233a9 100755 --- a/src/python/example/rips_complex_diagram_persistence_from_off_file_example.py +++ b/src/python/example/rips_complex_diagram_persistence_from_off_file_example.py @@ -3,7 +3,6 @@ import argparse import errno import os -import matplotlib.pyplot as plot import gudhi """ This file is part of the Gudhi Library - https://gudhi.inria.fr/ - @@ -70,6 +69,7 @@ with open(args.file, "r") as f: print(simplex_tree.betti_numbers()) if args.no_diagram == False: + import matplotlib.pyplot as plot gudhi.plot_persistence_diagram(diag, band=args.band) plot.show() else: diff --git a/src/python/example/rips_complex_edge_collapse_example.py b/src/python/example/rips_complex_edge_collapse_example.py new file mode 100755 index 00000000..b26eb9fc --- /dev/null +++ b/src/python/example/rips_complex_edge_collapse_example.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python + +import gudhi +import matplotlib.pyplot as plt +import time + +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + Author(s): Vincent Rouvreau + + Copyright (C) 2016 Inria + + Modification(s): + - YYYY/MM Author: Description of the modification +""" + +__author__ = "Vincent Rouvreau" +__copyright__ = "Copyright (C) 2020 Inria" +__license__ = "MIT" + + +print("#####################################################################") +print("RipsComplex (only the one-skeleton) creation from tore3D_300.off file") + +off_file = gudhi.__root_source_dir__ + '/data/points/tore3D_300.off' +point_cloud = gudhi.read_points_from_off_file(off_file = off_file) +rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=12.0) +simplex_tree = rips_complex.create_simplex_tree(max_dimension=1) +print('1. Rips complex is of dimension ', simplex_tree.dimension(), ' - ', + simplex_tree.num_simplices(), ' simplices - ', + simplex_tree.num_vertices(), ' vertices.') + +# Expansion of this one-skeleton would require a lot of memory. Let's collapse it +start = time.process_time() +simplex_tree.collapse_edges() +print('2. Rips complex is of dimension ', simplex_tree.dimension(), ' - ', + simplex_tree.num_simplices(), ' simplices - ', + simplex_tree.num_vertices(), ' vertices.') +simplex_tree.expansion(3) +diag = simplex_tree.persistence() +print("Collapse, expansion and persistence computation took ", time.process_time() - start, " sec.") + +# Use subplots to display diagram and density side by side +fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12, 5)) +gudhi.plot_persistence_diagram(diag, axes=axes[0]) +axes[0].set_title("Persistence after 1 collapse") + +# Collapse can be performed several times. Let's collapse it 3 times +start = time.process_time() +simplex_tree.collapse_edges(nb_iterations = 3) +print('3. Rips complex is of dimension ', simplex_tree.dimension(), ' - ', + simplex_tree.num_simplices(), ' simplices - ', + simplex_tree.num_vertices(), ' vertices.') +simplex_tree.expansion(3) +diag = simplex_tree.persistence() +print("Collapse, expansion and persistence computation took ", time.process_time() - start, " sec.") + +gudhi.plot_persistence_diagram(diag, axes=axes[1]) +axes[1].set_title("Persistence after 3 more collapses") + +# Plot the 2 persistence diagrams side to side to check the persistence is the same +plt.show()
\ No newline at end of file diff --git a/src/python/example/tangential_complex_plain_homology_from_off_file_example.py b/src/python/example/tangential_complex_plain_homology_from_off_file_example.py index 85bade4a..a4b4e9f5 100755 --- a/src/python/example/tangential_complex_plain_homology_from_off_file_example.py +++ b/src/python/example/tangential_complex_plain_homology_from_off_file_example.py @@ -3,7 +3,6 @@ import argparse import errno import os -import matplotlib.pyplot as plot import gudhi """ This file is part of the Gudhi Library - https://gudhi.inria.fr/ - @@ -62,6 +61,7 @@ with open(args.file, "r") as f: print(st.betti_numbers()) if args.no_diagram == False: + import matplotlib.pyplot as plot gudhi.plot_persistence_diagram(diag, band=args.band) plot.show() else: diff --git a/src/python/gudhi/__init__.py.in b/src/python/gudhi/__init__.py.in index 79e12fbc..3043201a 100644 --- a/src/python/gudhi/__init__.py.in +++ b/src/python/gudhi/__init__.py.in @@ -23,6 +23,10 @@ __all__ = [@GUDHI_PYTHON_MODULES@ @GUDHI_PYTHON_MODULES_EXTRA@] __available_modules = '' __missing_modules = '' +# For unitary tests purpose +# could use "if 'collapse_edges' in gudhi.__all__" when collapse edges will have a python module +__GUDHI_USE_EIGEN3 = @GUDHI_USE_EIGEN3@ + # Try to import * from gudhi.__module_name for default modules. # Extra modules require an explicit import by the user (mostly because of # unusual dependencies, but also to avoid cluttering namespace gudhi and diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx index ea128743..a4888914 100644 --- a/src/python/gudhi/alpha_complex.pyx +++ b/src/python/gudhi/alpha_complex.pyx @@ -16,7 +16,7 @@ from libcpp.utility cimport pair from libcpp.string cimport string from libcpp cimport bool from libc.stdint cimport intptr_t -import os +import warnings from gudhi.simplex_tree cimport * from gudhi.simplex_tree import SimplexTree @@ -28,66 +28,72 @@ __license__ = "GPL v3" cdef extern from "Alpha_complex_interface.h" namespace "Gudhi": cdef cppclass Alpha_complex_interface "Gudhi::alpha_complex::Alpha_complex_interface": - Alpha_complex_interface(vector[vector[double]] points, bool fast_version, bool exact_version) nogil except + + Alpha_complex_interface(vector[vector[double]] points, vector[double] weights, bool fast_version, bool exact_version) nogil except + vector[double] get_point(int vertex) nogil except + void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool default_filtration_value) nogil except + # AlphaComplex python interface cdef class AlphaComplex: - """AlphaComplex is a simplicial complex constructed from the finite cells - of a Delaunay Triangulation. + """AlphaComplex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation. - The filtration value of each simplex is computed as the square of the - circumradius of the simplex if the circumsphere is empty (the simplex is - then said to be Gabriel), and as the minimum of the filtration values of - the codimension 1 cofaces that make it not Gabriel otherwise. + The filtration value of each simplex is computed as the square of the circumradius of the simplex if the + circumsphere is empty (the simplex is then said to be Gabriel), and as the minimum of the filtration values of the + codimension 1 cofaces that make it not Gabriel otherwise. - All simplices that have a filtration value strictly greater than a given - alpha squared value are not inserted into the complex. + All simplices that have a filtration value strictly greater than a given alpha squared value are not inserted into + the complex. .. note:: - When Alpha_complex is constructed with an infinite value of alpha, the - complex is a Delaunay complex. - + When Alpha_complex is constructed with an infinite value of alpha, the complex is a Delaunay complex. """ cdef Alpha_complex_interface * this_ptr # Fake constructor that does nothing but documenting the constructor - def __init__(self, points=None, off_file='', precision='safe'): + def __init__(self, points=[], off_file='', weights=None, precision='safe'): """AlphaComplex constructor. :param points: A list of points in d-Dimension. - :type points: list of list of double - - Or + :type points: Iterable[Iterable[float]] - :param off_file: An OFF file style name. + :param off_file: **[deprecated]** An `OFF file style <fileformats.html#off-file-format>`_ name. + If an `off_file` is given with `points` as arguments, only points from the file are taken into account. :type off_file: string + :param weights: A list of weights. If set, the number of weights must correspond to the number of points. + :type weights: Iterable[float] + :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. :type precision: string + + :raises FileNotFoundError: **[deprecated]** If `off_file` is set but not found. + :raises ValueError: In case of inconsistency between the number of points and weights. """ # The real cython constructor - def __cinit__(self, points = None, off_file = '', precision = 'safe'): + def __cinit__(self, points = [], off_file = '', weights=None, precision = 'safe'): assert precision in ['fast', 'safe', 'exact'], "Alpha complex precision can only be 'fast', 'safe' or 'exact'" cdef bool fast = precision == 'fast' cdef bool exact = precision == 'exact' - cdef vector[vector[double]] pts if off_file: - if os.path.isfile(off_file): - points = read_points_from_off_file(off_file = off_file) - else: - print("file " + off_file + " not found.") - if points is None: - # Empty Alpha construction - points=[] + warnings.warn("off_file is a deprecated parameter, please consider using gudhi.read_points_from_off_file", + DeprecationWarning) + points = read_points_from_off_file(off_file = off_file) + + # weights are set but is inconsistent with the number of points + if weights != None and len(weights) != len(points): + raise ValueError("Inconsistency between the number of points and weights") + + # need to copy the points to use them without the gil + cdef vector[vector[double]] pts + cdef vector[double] wgts pts = points + if weights != None: + wgts = weights with nogil: - self.this_ptr = new Alpha_complex_interface(pts, fast, exact) + self.this_ptr = new Alpha_complex_interface(pts, wgts, fast, exact) def __dealloc__(self): if self.this_ptr != NULL: 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/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py index c6766c70..848dc03e 100644 --- a/src/python/gudhi/persistence_graphical_tools.py +++ b/src/python/gudhi/persistence_graphical_tools.py @@ -20,6 +20,7 @@ __author__ = "Vincent Rouvreau, Bertrand Michel, Theo Lacombe" __copyright__ = "Copyright (C) 2016 Inria" __license__ = "MIT" +_gudhi_matplotlib_use_tex = True def __min_birth_max_death(persistence, band=0.0): """This function returns (min_birth, max_death) from the persistence. @@ -117,10 +118,13 @@ def plot_persistence_barcode( try: import matplotlib.pyplot as plt import matplotlib.patches as mpatches - if _matplotlib_can_use_tex(): - from matplotlib import rc + from matplotlib import rc + if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex(): plt.rc('text', usetex=True) plt.rc('font', family='serif') + else: + plt.rc('text', usetex=False) + plt.rc('font', family='DejaVu Sans') if persistence_file != "": if path.isfile(persistence_file): @@ -263,10 +267,13 @@ def plot_persistence_diagram( try: import matplotlib.pyplot as plt import matplotlib.patches as mpatches - if _matplotlib_can_use_tex(): - from matplotlib import rc + from matplotlib import rc + if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex(): plt.rc('text', usetex=True) plt.rc('font', family='serif') + else: + plt.rc('text', usetex=False) + plt.rc('font', family='DejaVu Sans') if persistence_file != "": if path.isfile(persistence_file): @@ -436,10 +443,13 @@ def plot_persistence_density( import matplotlib.pyplot as plt import matplotlib.patches as mpatches from scipy.stats import kde - if _matplotlib_can_use_tex(): - from matplotlib import rc + from matplotlib import rc + if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex(): plt.rc('text', usetex=True) plt.rc('font', family='serif') + else: + plt.rc('text', usetex=False) + plt.rc('font', family='DejaVu Sans') if persistence_file != "": if dimension is None: diff --git a/src/python/gudhi/point_cloud/knn.py b/src/python/gudhi/point_cloud/knn.py index 4652fe80..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 @@ -46,7 +47,7 @@ class KNearestNeighbors: sort_results (bool): if True, then distances and indices of each point are sorted on return, so that the first column contains the closest points. Otherwise, neighbors are returned in an arbitrary order. Defaults to True. - enable_autodiff (bool): if the input is a torch.tensor, jax.numpy.ndarray or tensorflow.Tensor, this + enable_autodiff (bool): if the input is a torch.tensor or tensorflow.Tensor, this instructs the function to compute distances in a way that works with automatic differentiation. This is experimental, not supported for all metrics, and requires the package EagerPy. Defaults to False. @@ -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 46fee086..f8078d03 100644 --- a/src/python/gudhi/representations/vector_methods.py +++ b/src/python/gudhi/representations/vector_methods.py @@ -1,16 +1,20 @@ # 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): Mathieu Carrière +# Author(s): Mathieu Carrière, Martin Royer, Gard Spreemann # -# Copyright (C) 2018-2019 Inria +# Copyright (C) 2018-2020 Inria # # Modification(s): -# - YYYY/MM Author: Description of the modification +# - 2020/06 Martin: ATOL integration +# - 2020/12 Gard: A more flexible Betti curve class capable of computing exact curves. +# - 2021/11 Vincent Rouvreau: factorize _automatic_sample_range import numpy as np from sklearn.base import BaseEstimator, TransformerMixin +from sklearn.exceptions import NotFittedError from sklearn.preprocessing import MinMaxScaler, MaxAbsScaler from sklearn.neighbors import DistanceMetric +from sklearn.metrics import pairwise from .preprocessing import DiagramScaler, BirthPersistenceTransform @@ -44,10 +48,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): @@ -93,6 +101,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. @@ -118,10 +148,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): @@ -217,10 +244,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): @@ -284,77 +308,162 @@ class Silhouette(BaseEstimator, TransformerMixin): """ return self.fit_transform([diag])[0,:] + class BettiCurve(BaseEstimator, TransformerMixin): """ - This is a class for computing Betti curves from a list of persistence diagrams. A Betti curve is a 1D piecewise-constant function obtained from the rank function. It is sampled evenly on a given range and the vector of samples is returned. See https://www.researchgate.net/publication/316604237_Time_Series_Classification_via_Topological_Data_Analysis for more details. + Compute Betti curves from persistence diagrams. There are several modes of operation: with a given resolution (with or without a sample_range), with a predefined grid, and with none of the previous. With a predefined grid, the class computes the Betti numbers at those grid points. Without a predefined grid, if the resolution is set to None, it can be fit to a list of persistence diagrams and produce a grid that consists of (at least) the filtration values at which at least one of those persistence diagrams changes Betti numbers, and then compute the Betti numbers at those grid points. In the latter mode, the exact Betti curve is computed for the entire real line. Otherwise, if the resolution is given, the Betti curve is obtained by sampling evenly using either the given sample_range or based on the persistence diagrams. """ - def __init__(self, resolution=100, sample_range=[np.nan, np.nan]): + + def __init__(self, resolution=100, sample_range=[np.nan, np.nan], predefined_grid=None): """ Constructor for the BettiCurve class. Parameters: resolution (int): number of sample for the piecewise-constant function (default 100). sample_range ([double, double]): minimum and maximum of the piecewise-constant function domain, of the form [x_min, x_max] (default [numpy.nan, numpy.nan]). It is the interval on which samples will be drawn evenly. If one of the values is numpy.nan, it can be computed from the persistence diagrams with the fit() method. + predefined_grid (1d array or None, default=None): Predefined filtration grid points at which to compute the Betti curves. Must be strictly ordered. Infinities are ok. If None (default), and resolution is given, the grid will be uniform from x_min to x_max in 'resolution' steps, otherwise a grid will be computed that captures all changes in Betti numbers in the provided data. + + Attributes: + grid_ (1d array): The grid on which the Betti numbers are computed. If predefined_grid was specified, `grid_` will always be that grid, independently of data. If not, the grid is fitted to capture all filtration values at which the Betti numbers change. + + Examples + -------- + If pd is a persistence diagram and xs is a nonempty grid of finite values such that xs[0] >= pd.min(), then the results of: + + >>> bc = BettiCurve(predefined_grid=xs) # doctest: +SKIP + >>> result = bc(pd) # doctest: +SKIP + + and + + >>> from scipy.interpolate import interp1d # doctest: +SKIP + >>> bc = BettiCurve(resolution=None, predefined_grid=None) # doctest: +SKIP + >>> bettis = bc.fit_transform([pd]) # doctest: +SKIP + >>> interp = interp1d(bc.grid_, bettis[0, :], kind="previous", fill_value="extrapolate") # doctest: +SKIP + >>> result = np.array(interp(xs), dtype=int) # doctest: +SKIP + + are the same. """ - self.resolution, self.sample_range = resolution, sample_range - def fit(self, X, y=None): + if (predefined_grid is not None) and (not isinstance(predefined_grid, np.ndarray)): + raise ValueError("Expected predefined_grid as array or None.") + + self.predefined_grid = predefined_grid + self.resolution = resolution + self.sample_range = sample_range + + def is_fitted(self): + return hasattr(self, "grid_") + + def fit(self, X, y = None): """ - Fit the BettiCurve class on a list of persistence diagrams: if any of the values in **sample_range** is numpy.nan, replace it with the corresponding value computed on the given list of persistence diagrams. + Fit the BettiCurve class on a list of persistence diagrams: if any of the values in **sample_range** is numpy.nan, replace it with the corresponding value computed on the given list of persistence diagrams. When no predefined grid is provided and resolution set to None, compute a filtration grid that captures all changes in Betti numbers for all the given persistence diagrams. Parameters: - X (list of n x 2 numpy arrays): input persistence diagrams. - y (n x 1 array): persistence diagram labels (unused). + X (list of 2d arrays): Persistence diagrams. + y (None): Ignored. """ - 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)) + + if self.predefined_grid is None: + if self.resolution is None: # Flexible/exact version + events = np.unique(np.concatenate([pd.flatten() for pd in X] + [[-np.inf]], axis=0)) + self.grid_ = np.array(events) + else: + self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y) + self.grid_ = np.linspace(self.sample_range[0], self.sample_range[1], self.resolution) + else: + self.grid_ = self.predefined_grid # Get the predefined grid from user + return self def transform(self, X): """ - Compute the Betti curve for each persistence diagram individually and concatenate the results. + Compute Betti curves. Parameters: - X (list of n x 2 numpy arrays): input persistence diagrams. - + X (list of 2d arrays): Persistence diagrams. + Returns: - numpy array with shape (number of diagrams) x (**resolution**): output Betti curves. + `len(X).len(self.grid_)` array of ints: Betti numbers of the given persistence diagrams at the grid points given in `self.grid_` """ - num_diag, Xfit = len(X), [] - x_values = np.linspace(self.sample_range[0], self.sample_range[1], self.resolution) - step_x = x_values[1] - x_values[0] - for i in range(num_diag): + if not self.is_fitted(): + raise NotFittedError("Not fitted.") - diagram, num_pts_in_diag = X[i], X[i].shape[0] + if not X: + X = [np.zeros((0, 2))] + + N = len(X) - bc = np.zeros(self.resolution) - for j in range(num_pts_in_diag): - [px,py] = diagram[j,:2] - min_idx = np.clip(np.ceil((px - self.sample_range[0]) / step_x).astype(int), 0, self.resolution) - 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): - bc[k] += 1 + events = np.concatenate([pd.flatten(order="F") for pd in X], axis=0) + sorting = np.argsort(events) + offsets = np.zeros(1 + N, dtype=int) + for i in range(0, N): + offsets[i+1] = offsets[i] + 2*X[i].shape[0] + starts = offsets[0:N] + ends = offsets[1:N + 1] - 1 - Xfit.append(np.reshape(bc,[1,-1])) + bettis = [[0] for i in range(0, N)] - Xfit = np.concatenate(Xfit, 0) + i = 0 + for x in self.grid_: + while i < len(sorting) and events[sorting[i]] <= x: + j = np.searchsorted(ends, sorting[i]) + delta = 1 if sorting[i] - starts[j] < len(X[j]) else -1 + bettis[j][-1] += delta + i += 1 + for k in range(0, N): + bettis[k].append(bettis[k][-1]) - return Xfit + return np.array(bettis, dtype=int)[:, 0:-1] - def __call__(self, diag): + def fit_transform(self, X): + """ + The result is the same as fit(X) followed by transform(X), but potentially faster. """ - Apply BettiCurve on a single persistence diagram and outputs the result. - Parameters: - diag (n x 2 numpy array): input persistence diagram. + if self.predefined_grid is None and self.resolution is None: + if not X: + X = [np.zeros((0, 2))] - Returns: - numpy array with shape (**resolution**): output Betti curve. + N = len(X) + + events = np.concatenate([pd.flatten(order="F") for pd in X], axis=0) + sorting = np.argsort(events) + offsets = np.zeros(1 + N, dtype=int) + for i in range(0, N): + offsets[i+1] = offsets[i] + 2*X[i].shape[0] + starts = offsets[0:N] + ends = offsets[1:N + 1] - 1 + + xs = [-np.inf] + bettis = [[0] for i in range(0, N)] + + for i in sorting: + j = np.searchsorted(ends, i) + delta = 1 if i - starts[j] < len(X[j]) else -1 + if events[i] == xs[-1]: + bettis[j][-1] += delta + else: + xs.append(events[i]) + for k in range(0, j): + bettis[k].append(bettis[k][-1]) + bettis[j].append(bettis[j][-1] + delta) + for k in range(j+1, N): + bettis[k].append(bettis[k][-1]) + + self.grid_ = np.array(xs) + return np.array(bettis, dtype=int) + + else: + return self.fit(X).transform(X) + + def __call__(self, diag): """ - return self.fit_transform([diag])[0,:] + Shorthand for transform on a single persistence diagram. + """ + return self.fit_transform([diag])[0, :] + + class Entropy(BaseEstimator, TransformerMixin): """ @@ -380,10 +489,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): @@ -402,9 +508,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])) ) @@ -418,12 +528,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): @@ -484,7 +593,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] @@ -574,3 +689,141 @@ class ComplexPolynomial(BaseEstimator, TransformerMixin): numpy array with shape (**threshold**): output complex vector of coefficients. """ return self.fit_transform([diag])[0,:] + +def _lapl_contrast(measure, centers, inertias): + """contrast function for vectorising `measure` in ATOL""" + return np.exp(-pairwise.pairwise_distances(measure, Y=centers) / inertias) + +def _gaus_contrast(measure, centers, inertias): + """contrast function for vectorising `measure` in ATOL""" + return np.exp(-pairwise.pairwise_distances(measure, Y=centers, squared=True) / inertias**2) + +def _indicator_contrast(diags, centers, inertias): + """contrast function for vectorising `measure` in ATOL""" + robe_curve = np.clip(2-pairwise.pairwise_distances(diags, Y=centers)/inertias, 0, 1) + return robe_curve + +def _cloud_weighting(measure): + """automatic uniform weighting with mass 1 for `measure` in ATOL""" + return np.ones(shape=measure.shape[0]) + +def _iidproba_weighting(measure): + """automatic uniform weighting with mass 1/N for `measure` in ATOL""" + return np.ones(shape=measure.shape[0]) / measure.shape[0] + +class Atol(BaseEstimator, TransformerMixin): + """ + This class allows to vectorise measures (e.g. point clouds, persistence diagrams, etc) after a quantisation step. + + ATOL paper: :cite:`royer2019atol` + + Example + -------- + >>> from sklearn.cluster import KMeans + >>> from gudhi.representations.vector_methods import Atol + >>> import numpy as np + >>> a = np.array([[1, 2, 4], [1, 4, 0], [1, 0, 4]]) + >>> b = np.array([[4, 2, 0], [4, 4, 0], [4, 0, 2]]) + >>> c = np.array([[3, 2, -1], [1, 2, -1]]) + >>> atol_vectoriser = Atol(quantiser=KMeans(n_clusters=2, random_state=202006)) + >>> atol_vectoriser.fit(X=[a, b, c]).centers # doctest: +SKIP + >>> # array([[ 2. , 0.66666667, 3.33333333], + >>> # [ 2.6 , 2.8 , -0.4 ]]) + >>> atol_vectoriser(a) + >>> # array([1.18168665, 0.42375966]) # doctest: +SKIP + >>> atol_vectoriser(c) + >>> # array([0.02062512, 1.25157463]) # doctest: +SKIP + >>> atol_vectoriser.transform(X=[a, b, c]) # doctest: +SKIP + >>> # array([[1.18168665, 0.42375966], + >>> # [0.29861028, 1.06330156], + >>> # [0.02062512, 1.25157463]]) + """ + # Note the example above must be up to date with the one in tests called test_atol_doc + def __init__(self, quantiser, weighting_method="cloud", contrast="gaussian"): + """ + Constructor for the Atol measure vectorisation class. + + Parameters: + quantiser (Object): Object with `fit` (sklearn API consistent) and `cluster_centers` and `n_clusters` + attributes, e.g. sklearn.cluster.KMeans. It will be fitted when the Atol object function `fit` is called. + weighting_method (string): constant generic function for weighting the measure points + choose from {"cloud", "iidproba"} + (default: constant function, i.e. the measure is seen as a point cloud by default). + This will have no impact if weights are provided along with measures all the way: `fit` and `transform`. + contrast (string): constant function for evaluating proximity of a measure with respect to centers + choose from {"gaussian", "laplacian", "indicator"} + (default: gaussian contrast function, see page 3 in the ATOL paper). + """ + self.quantiser = quantiser + self.contrast = { + "gaussian": _gaus_contrast, + "laplacian": _lapl_contrast, + "indicator": _indicator_contrast, + }.get(contrast, _gaus_contrast) + self.weighting_method = { + "cloud" : _cloud_weighting, + "iidproba": _iidproba_weighting, + }.get(weighting_method, _cloud_weighting) + + def fit(self, X, y=None, sample_weight=None): + """ + Calibration step: fit centers to the sample measures and derive inertias between centers. + + Parameters: + X (list N x d numpy arrays): input measures in R^d from which to learn center locations and inertias + (measures can have different N). + y: Ignored, present for API consistency by convention. + sample_weight (list of numpy arrays): weights for each measure point in X, optional. + If None, the object's weighting_method will be used. + + Returns: + self + """ + if not hasattr(self.quantiser, 'fit'): + raise TypeError("quantiser %s has no `fit` attribute." % (self.quantiser)) + if sample_weight is None: + sample_weight = np.concatenate([self.weighting_method(measure) for measure in X]) + + measures_concat = np.concatenate(X) + self.quantiser.fit(X=measures_concat, sample_weight=sample_weight) + self.centers = self.quantiser.cluster_centers_ + if self.quantiser.n_clusters == 1: + dist_centers = pairwise.pairwise_distances(measures_concat) + np.fill_diagonal(dist_centers, 0) + self.inertias = np.array([np.max(dist_centers)/2]) + else: + dist_centers = pairwise.pairwise_distances(self.centers) + dist_centers[dist_centers == 0] = np.inf + self.inertias = np.min(dist_centers, axis=0)/2 + return self + + def __call__(self, measure, sample_weight=None): + """ + Apply measure vectorisation on a single measure. + + Parameters: + measure (n x d numpy array): input measure in R^d. + + Returns: + numpy array in R^self.quantiser.n_clusters. + """ + if sample_weight is None: + sample_weight = self.weighting_method(measure) + return np.sum(sample_weight * self.contrast(measure, self.centers, self.inertias.T).T, axis=1) + + def transform(self, X, sample_weight=None): + """ + Apply measure vectorisation on a list of measures. + + Parameters: + X (list N x d numpy arrays): input measures in R^d from which to learn center locations and inertias + (measures can have different N). + sample_weight (list of numpy arrays): weights for each measure point in X, optional. + If None, the object's weighting_method will be used. + + Returns: + numpy array with shape (number of measures) x (self.quantiser.n_clusters). + """ + if sample_weight is None: + sample_weight = [self.weighting_method(measure) for measure in X] + return np.stack([self(measure, sample_weight=weight) for measure, weight in zip(X, sample_weight)]) diff --git a/src/python/gudhi/rips_complex.pyx b/src/python/gudhi/rips_complex.pyx index 72e82c79..c3470292 100644 --- a/src/python/gudhi/rips_complex.pyx +++ b/src/python/gudhi/rips_complex.pyx @@ -49,13 +49,13 @@ cdef class RipsComplex: :type max_edge_length: float :param points: A list of points in d-Dimension. - :type points: list of list of double + :type points: list of list of float Or :param distance_matrix: A distance matrix (full square or lower triangular). - :type points: list of list of double + :type points: list of list of float And in both cases @@ -89,10 +89,10 @@ cdef class RipsComplex: def create_simplex_tree(self, max_dimension=1): """ - :param max_dimension: graph expansion for rips until this given maximal + :param max_dimension: graph expansion for Rips until this given maximal dimension. :type max_dimension: int - :returns: A simplex tree created from the Delaunay Triangulation. + :returns: A simplex tree encoding the Vietoris–Rips filtration. :rtype: SimplexTree """ stree = SimplexTree() diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 3bd9b080..175f970e 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -36,9 +36,16 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": Simplex_tree_skeleton_iterator operator++() nogil bint operator!=(Simplex_tree_skeleton_iterator) nogil + cdef cppclass Simplex_tree_boundary_iterator "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>::Boundary_simplex_iterator": + Simplex_tree_boundary_iterator() nogil + Simplex_tree_simplex_handle& operator*() nogil + Simplex_tree_boundary_iterator operator++() nogil + bint operator!=(Simplex_tree_boundary_iterator) nogil + 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 + Simplex_tree_interface_full_featured(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 @@ -58,6 +65,9 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": bool make_filtration_non_decreasing() nogil void compute_extended_filtration() nogil vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]] dgm, double min_persistence) nogil + Simplex_tree_interface_full_featured* collapse_edges(int nb_collapse_iteration) nogil except + + void reset_filtration(double filtration, int dimension) nogil + bint operator==(Simplex_tree_interface_full_featured) nogil # Iterators over Simplex tree pair[vector[int], double] get_simplex_and_filtration(Simplex_tree_simplex_handle f_simplex) nogil Simplex_tree_simplices_iterator get_simplices_iterator_begin() nogil @@ -66,11 +76,12 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": vector[Simplex_tree_simplex_handle].const_iterator get_filtration_iterator_end() nogil Simplex_tree_skeleton_iterator get_skeleton_iterator_begin(int dimension) nogil Simplex_tree_skeleton_iterator get_skeleton_iterator_end(int dimension) nogil + pair[Simplex_tree_boundary_iterator, Simplex_tree_boundary_iterator] get_boundary_iterators(vector[int] simplex) nogil except + 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 42f1e635..3eebfff7 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" @@ -39,13 +38,29 @@ cdef class SimplexTree: cdef Simplex_tree_persistence_interface * pcohptr # Fake constructor that does nothing but documenting the constructor - def __init__(self): + def __init__(self, other = None): """SimplexTree constructor. + + :param other: If `other` is `None` (default value), an empty `SimplexTree` is created. + If `other` is a `SimplexTree`, the `SimplexTree` is constructed from a deep copy of `other`. + :type other: SimplexTree (Optional) + :returns: An empty or a copy simplex tree. + :rtype: SimplexTree + + :raises TypeError: In case `other` is neither `None`, nor a `SimplexTree`. + :note: If the `SimplexTree` is a copy, the persistence information is not copied. If you need it in the clone, + you have to call :func:`compute_persistence` on it even if you had already computed it in the original. """ # The real cython constructor - def __cinit__(self): - self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured()) + def __cinit__(self, other = None): + if other: + if isinstance(other, SimplexTree): + self.thisptr = _get_copy_intptr(other) + else: + raise TypeError("`other` argument requires to be of type `SimplexTree`, or `None`.") + else: + self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured()) def __dealloc__(self): cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr() @@ -64,12 +79,27 @@ cdef class SimplexTree: """ return self.pcohptr != NULL + def copy(self): + """ + :returns: A simplex tree that is a deep copy of itself. + :rtype: SimplexTree + + :note: The persistence information is not copied. If you need it in the clone, you have to call + :func:`compute_persistence` on it even if you had already computed it in the original. + """ + stree = SimplexTree() + stree.thisptr = _get_copy_intptr(self) + return stree + + def __deepcopy__(self): + return self.copy() + def filtration(self, simplex): """This function returns the filtration value for a given N-simplex in this simplicial complex, or +infinity if it is not in the complex. :param simplex: The N-simplex, represented by a list of vertex. - :type simplex: list of int. + :type simplex: list of int :returns: The simplicial complex filtration value. :rtype: float """ @@ -80,7 +110,7 @@ cdef class SimplexTree: given N-simplex. :param simplex: The N-simplex, represented by a list of vertex. - :type simplex: list of int. + :type simplex: list of int :param filtration: The new filtration value. :type filtration: float @@ -153,7 +183,7 @@ cdef class SimplexTree: """This function sets the dimension of the simplicial complex. :param dimension: The new dimension value. - :type dimension: int. + :type dimension: int .. note:: @@ -172,7 +202,7 @@ cdef class SimplexTree: complex or not. :param simplex: The N-simplex to find, represented by a list of vertex. - :type simplex: list of int. + :type simplex: list of int :returns: true if the simplex was found, false otherwise. :rtype: bool """ @@ -186,9 +216,9 @@ cdef class SimplexTree: :param simplex: The N-simplex to insert, represented by a list of vertex. - :type simplex: list of int. + :type simplex: list of int :param filtration: The filtration value of the simplex. - :type filtration: float. + :type filtration: float :returns: true if the simplex was not yet in the complex, false otherwise (whatever its original filtration value). :rtype: bool @@ -196,7 +226,7 @@ cdef class SimplexTree: return self.get_ptr().insert(simplex, <double>filtration) @staticmethod - def create_from_array(filtrations, double max_filtration=numpy.inf): + def create_from_array(filtrations, double max_filtration=np.inf): """Creates a new, empty complex and inserts vertices and edges. The vertices are numbered from 0 to n-1, and the filtration values are encoded in the array, with the diagonal representing the vertices. It is the caller's responsibility to ensure that this defines a filtration, which can be achieved with either:: @@ -216,7 +246,7 @@ cdef class SimplexTree: :rtype: SimplexTree """ # TODO: document which half of the matrix is actually read? - filtrations = numpy.asanyarray(filtrations, dtype=float) + filtrations = np.asanyarray(filtrations, dtype=float) cdef double[:,:] F = filtrations ret = SimplexTree() cdef int n = F.shape[0] @@ -272,7 +302,7 @@ cdef class SimplexTree: """This function returns a generator with the (simplices of the) skeleton of a maximum given dimension. :param dimension: The skeleton dimension value. - :type dimension: int. + :type dimension: int :returns: The (simplices of the) skeleton of a maximum dimension. :rtype: generator with tuples(simplex, filtration) """ @@ -287,7 +317,7 @@ cdef class SimplexTree: """This function returns the star of a given N-simplex. :param simplex: The N-simplex, represented by a list of vertex. - :type simplex: list of int. + :type simplex: list of int :returns: The (simplices of the) star of a simplex. :rtype: list of tuples(simplex, filtration) """ @@ -309,10 +339,10 @@ cdef class SimplexTree: given codimension. :param simplex: The N-simplex, represented by a list of vertex. - :type simplex: list of int. + :type simplex: list of int :param codimension: The codimension. If codimension = 0, all cofaces are returned (equivalent of get_star function) - :type codimension: int. + :type codimension: int :returns: The (simplices of the) cofaces of a simplex :rtype: list of tuples(simplex, filtration) """ @@ -329,12 +359,28 @@ cdef class SimplexTree: ct.append((v, filtered_simplex.second)) return ct + def get_boundaries(self, simplex): + """This function returns a generator with the boundaries of a given N-simplex. + If you do not need the filtration values, the boundary can also be obtained as + :code:`itertools.combinations(simplex,len(simplex)-1)`. + + :param simplex: The N-simplex, represented by a list of vertex. + :type simplex: list of int. + :returns: The (simplices of the) boundary of a simplex + :rtype: generator with tuples(simplex, filtration) + """ + cdef pair[Simplex_tree_boundary_iterator, Simplex_tree_boundary_iterator] it = self.get_ptr().get_boundary_iterators(simplex) + + while it.first != it.second: + yield self.get_ptr().get_simplex_and_filtration(dereference(it.first)) + preincrement(it.first) + def remove_maximal_simplex(self, simplex): """This function removes a given maximal N-simplex from the simplicial complex. :param simplex: The N-simplex, represented by a list of vertex. - :type simplex: list of int. + :type simplex: list of int .. note:: @@ -352,7 +398,7 @@ cdef class SimplexTree: """Prune above filtration value given as parameter. :param filtration: Maximum threshold value. - :type filtration: float. + :type filtration: float :returns: The filtration modification information. :rtype: bool @@ -372,7 +418,7 @@ cdef class SimplexTree: return self.get_ptr().prune_above_filtration(filtration) def expansion(self, max_dim): - """Expands the Simplex_tree containing only its one skeleton + """Expands the simplex tree containing only its one skeleton until dimension max_dim. The expanded simplicial complex until dimension :math:`d` @@ -382,11 +428,11 @@ cdef class SimplexTree: The filtration value assigned to a simplex is the maximal filtration value of one of its edges. - The Simplex_tree must contain no simplex of dimension bigger than + The simplex tree must contain no simplex of dimension bigger than 1 when calling the method. :param max_dim: The maximal dimension. - :type max_dim: int. + :type max_dim: int """ cdef int maxdim = max_dim with nogil: @@ -402,38 +448,54 @@ cdef class SimplexTree: """ return self.get_ptr().make_filtration_non_decreasing() + def reset_filtration(self, filtration, min_dim = 0): + """This function resets the filtration value of all the simplices of dimension at least min_dim. Resets all the + simplex tree when `min_dim = 0`. + `reset_filtration` may break the filtration property with `min_dim > 0`, and it is the user's responsibility to + make it a valid filtration (using a large enough `filt_value`, or calling `make_filtration_non_decreasing` + afterwards for instance). + + :param filtration: New threshold value. + :type filtration: float. + :param min_dim: The minimal dimension. Default value is 0. + :type min_dim: int. + """ + self.get_ptr().reset_filtration(filtration, min_dim) + def extend_filtration(self): - """ Extend filtration for computing extended persistence. This function only uses the - filtration values at the 0-dimensional simplices, and computes the extended persistence - diagram induced by the lower-star filtration computed with these values. + """ Extend filtration for computing extended persistence. This function only uses the filtration values at the + 0-dimensional simplices, and computes the extended persistence diagram induced by the lower-star filtration + computed with these values. .. note:: - Note that after calling this function, the filtration - values are actually modified within the Simplex_tree. - The function :func:`extended_persistence` - retrieves the original values. + Note that after calling this function, the filtration values are actually modified within the simplex tree. + The function :func:`extended_persistence` retrieves the original values. .. note:: - Note that this code creates an extra vertex internally, so you should make sure that - the Simplex_tree does not contain a vertex with the largest possible value (i.e., 4294967295). + Note that this code creates an extra vertex internally, so you should make sure that the simplex tree does + not contain a vertex with the largest possible value (i.e., 4294967295). + + This `notebook <https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-extended-persistence.ipynb>`_ + explains how to compute an extension of persistence called extended persistence. """ self.get_ptr().compute_extended_filtration() def extended_persistence(self, homology_coeff_field=11, min_persistence=0): - """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. - :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). Default value is - 0.0. - Sets min_persistence to -1.0 to see all values. - :type min_persistence: float. - :returns: A list of four persistence diagrams in the format described in :func:`persistence`. The first one is Ordinary, the second one is Relative, the third one is Extended+ and the fourth one is Extended-. See https://link.springer.com/article/10.1007/s10208-008-9027-z and/or section 2.2 in https://link.springer.com/article/10.1007/s10208-017-9370-z for a description of these subtypes. + """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. 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). + Default value is 0.0. Sets min_persistence to -1.0 to see all values. + :type min_persistence: float + :returns: A list of four persistence diagrams in the format described in :func:`persistence`. The first one is + Ordinary, the second one is Relative, the third one is Extended+ and the fourth one is Extended-. + See https://link.springer.com/article/10.1007/s10208-008-9027-z and/or section 2.2 in + https://link.springer.com/article/10.1007/s10208-017-9370-z for a description of these subtypes. .. note:: @@ -444,6 +506,9 @@ cdef class SimplexTree: The coordinates of the persistence diagram points might be a little different than the original filtration values due to the internal transformation (scaling to [-2,-1]) that is performed on these values during the computation of extended persistence. + + This `notebook <https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-extended-persistence.ipynb>`_ + explains how to compute an extension of persistence called extended persistence. """ cdef vector[pair[int, pair[double, double]]] persistence_result if self.pcohptr != NULL: @@ -458,13 +523,13 @@ 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. - :type homology_coeff_field: int. + 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 0.0. Set min_persistence to -1.0 to see all values. - :type min_persistence: float. + :type min_persistence: float :param persistence_dim_max: If true, the persistent homology for the maximal dimension in the complex is computed. If false, it is ignored. Default is false. @@ -481,13 +546,13 @@ 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. - :type homology_coeff_field: int. + 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 0.0. Sets min_persistence to -1.0 to see all values. - :type min_persistence: float. + :type min_persistence: float :param persistence_dim_max: If true, the persistent homology for the maximal dimension in the complex is computed. If false, it is ignored. Default is false. @@ -522,10 +587,10 @@ cdef class SimplexTree: :param from_value: The persistence birth limit to be added in the numbers (persistent birth <= from_value). - :type from_value: float. + :type from_value: float :param to_value: The persistence death limit to be added in the numbers (persistent death > to_value). - :type to_value: float. + :type to_value: float :returns: The persistent Betti numbers ([B0, B1, ..., Bn]). :rtype: list of int @@ -542,7 +607,7 @@ cdef class SimplexTree: complex in a specific dimension. :param dimension: The specific dimension. - :type dimension: int. + :type dimension: int :returns: The persistence intervals. :rtype: numpy array of dimension 2 @@ -551,7 +616,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. @@ -571,7 +640,7 @@ cdef class SimplexTree: complex in a user given file name. :param persistence_file: Name of the file. - :type persistence_file: string. + :type persistence_file: string :note: intervals_in_dim function requires :func:`compute_persistence` @@ -592,8 +661,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): @@ -611,17 +680,50 @@ 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): + """Assuming the simplex tree is a 1-skeleton graph, this method collapse edges (simplices of higher dimension + are ignored) and resets the simplex tree from the remaining edges. + A good candidate is to build a simplex tree on top of a :class:`~gudhi.RipsComplex` of dimension 1 before + collapsing edges + (cf. :download:`rips_complex_edge_collapse_example.py <../example/rips_complex_edge_collapse_example.py>`). + For implementation details, please refer to :cite:`edgecollapsesocg2020`. + + :param nb_iterations: The number of edge collapse iterations to perform. Default is 1. + :type nb_iterations: int + + :note: collapse_edges method requires `Eigen <installation.html#eigen>`_ >= 3.1.0 and an exception is thrown + if this method is not available. + """ + # Backup old pointer + cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr() + cdef int nb_iter = nb_iterations + with nogil: + # New pointer is a new collapsed simplex tree + self.thisptr = <intptr_t>(ptr.collapse_edges(nb_iter)) + # Delete old pointer + del ptr + + def __eq__(self, other:SimplexTree): + """Test for structural equality + :returns: True if the 2 simplex trees are equal, False otherwise. + :rtype: bool + """ + return dereference(self.get_ptr()) == dereference(other.get_ptr()) + +cdef intptr_t _get_copy_intptr(SimplexTree stree) nogil: + return <intptr_t>(new Simplex_tree_interface_full_featured(dereference(stree.get_ptr()))) diff --git a/src/python/gudhi/subsampling.pyx b/src/python/gudhi/subsampling.pyx index f77c6f75..46f32335 100644 --- a/src/python/gudhi/subsampling.pyx +++ b/src/python/gudhi/subsampling.pyx @@ -33,7 +33,7 @@ def choose_n_farthest_points(points=None, off_file='', nb_points=0, starting_poi The iteration starts with the landmark `starting point`. :param points: The input point set. - :type points: Iterable[Iterable[float]]. + :type points: Iterable[Iterable[float]] Or @@ -42,14 +42,15 @@ def choose_n_farthest_points(points=None, off_file='', nb_points=0, starting_poi And in both cases - :param nb_points: Number of points of the subsample. - :type nb_points: unsigned. + :param nb_points: Number of points of the subsample (the subsample may be \ + smaller if there are fewer than nb_points distinct input points) + :type nb_points: int :param starting_point: The iteration starts with the landmark `starting \ - point`,which is the index of the point to start with. If not set, this \ + point`, which is the index of the point to start with. If not set, this \ index is chosen randomly. - :type starting_point: unsigned. + :type starting_point: int :returns: The subsample point set. - :rtype: List[List[float]]. + :rtype: List[List[float]] """ if off_file: if os.path.isfile(off_file): @@ -76,7 +77,7 @@ def pick_n_random_points(points=None, off_file='', nb_points=0): """Subsample a point set by picking random vertices. :param points: The input point set. - :type points: Iterable[Iterable[float]]. + :type points: Iterable[Iterable[float]] Or @@ -86,7 +87,7 @@ def pick_n_random_points(points=None, off_file='', nb_points=0): And in both cases :param nb_points: Number of points of the subsample. - :type nb_points: unsigned. + :type nb_points: int :returns: The subsample point set. :rtype: List[List[float]] """ @@ -104,10 +105,10 @@ def pick_n_random_points(points=None, off_file='', nb_points=0): def sparsify_point_set(points=None, off_file='', min_squared_dist=0.0): """Outputs a subset of the input points so that the squared distance - between any two points is greater than or equal to min_squared_dist. + between any two points is greater than min_squared_dist. :param points: The input point set. - :type points: Iterable[Iterable[float]]. + :type points: Iterable[Iterable[float]] Or @@ -118,7 +119,7 @@ def sparsify_point_set(points=None, off_file='', min_squared_dist=0.0): :param min_squared_dist: Minimum squared distance separating the output \ points. - :type min_squared_dist: float. + :type min_squared_dist: float :returns: The subsample point set. :rtype: List[List[float]] """ diff --git a/src/python/gudhi/wasserstein/wasserstein.py b/src/python/gudhi/wasserstein/wasserstein.py index b37d30bb..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, tensorflow.Tensor or jax.numpy.ndarray, 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) @@ -165,17 +339,17 @@ def wasserstein_distance(X, Y, matching=False, order=1., internal_p=np.inf, enab # empty arrays are not handled properly by the helpers, so we avoid calling them if len(pairs_X_Y): dists.append((Y_orig[pairs_X_Y[:, 1]] - X_orig[pairs_X_Y[:, 0]]).norms.lp(internal_p, axis=-1).norms.lp(order)) - if len(pairs_X_diag): + if len(pairs_X_diag[0]): dists.append(_perstot_autodiff(X_orig[pairs_X_diag], order, internal_p)) - if len(pairs_Y_diag): + if len(pairs_Y_diag[0]): dists.append(_perstot_autodiff(Y_orig[pairs_Y_diag], order, internal_p)) dists = [dist.reshape(1) for dist in dists] 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/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py index 0541572b..16f63c3d 100644 --- a/src/python/gudhi/weighted_rips_complex.py +++ b/src/python/gudhi/weighted_rips_complex.py @@ -12,9 +12,11 @@ from gudhi import SimplexTree class WeightedRipsComplex: """ Class to generate a weighted Rips complex from a distance matrix and weights on vertices, - in the way described in :cite:`dtmfiltrations`. + in the way described in :cite:`dtmfiltrations` with `p=1`. The filtration value of vertex `i` is `2*weights[i]`, + and the filtration value of edge `ij` is `distance_matrix[i][j]+weights[i]+weights[j]`, + or the maximum of the filtrations of its extremities, whichever is largest. Remark that all the filtration values are doubled compared to the definition in the paper - for the consistency with RipsComplex. + for consistency with RipsComplex. """ def __init__(self, distance_matrix, diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index d699ad9b..3d20aa8f 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -31,15 +31,34 @@ namespace Gudhi { namespace alpha_complex { -template <typename CgalPointType> -std::vector<double> pt_cgal_to_cython(CgalPointType const& point) { - std::vector<double> vd; - vd.reserve(point.dimension()); - for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) - vd.push_back(CGAL::to_double(*coord)); - return vd; -} +// template Functor that transforms a CGAL point to a vector of double as expected by cython +template<typename CgalPointType, bool Weighted> +struct Point_cgal_to_cython; + +// Specialized Unweighted Functor +template<typename CgalPointType> +struct Point_cgal_to_cython<CgalPointType, false> { + std::vector<double> operator()(CgalPointType const& point) const + { + std::vector<double> vd; + vd.reserve(point.dimension()); + for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) + vd.push_back(CGAL::to_double(*coord)); + return vd; + } +}; +// Specialized Weighted Functor +template<typename CgalPointType> +struct Point_cgal_to_cython<CgalPointType, true> { + std::vector<double> operator()(CgalPointType const& weighted_point) const + { + const auto& point = weighted_point.point(); + return Point_cgal_to_cython<decltype(point), false>()(point); + } +}; + +// Function that transforms a cython point (aka. a vector of double) to a CGAL point template <typename CgalPointType> static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) { return CgalPointType(vec.size(), vec.begin(), vec.end()); @@ -48,24 +67,38 @@ static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) { class Abstract_alpha_complex { public: virtual std::vector<double> get_point(int vh) = 0; + virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) = 0; + + virtual std::size_t num_vertices() const = 0; + + virtual ~Abstract_alpha_complex() = default; }; -class Exact_Alphacomplex_dD : public Abstract_alpha_complex { +template <bool Weighted = false> +class Exact_alpha_complex_dD final : public Abstract_alpha_complex { private: using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; - using Point = typename Kernel::Point_d; + using Bare_point = typename Kernel::Point_d; + using Point = std::conditional_t<Weighted, typename Kernel::Weighted_point_d, + typename Kernel::Point_d>; public: - Exact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version) + Exact_alpha_complex_dD(const std::vector<std::vector<double>>& points, bool exact_version) : exact_version_(exact_version), - alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>)) { + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>)) { + } + + Exact_alpha_complex_dD(const std::vector<std::vector<double>>& points, + const std::vector<double>& weights, bool exact_version) + : exact_version_(exact_version), + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>), weights) { } virtual std::vector<double> get_point(int vh) override { - Point const& point = alpha_complex_.get_point(vh); - return pt_cgal_to_cython(point); + // Can be a Weighted or a Bare point in function of Weighted + return Point_cgal_to_cython<Point, Weighted>()(alpha_complex_.get_point(vh)); } virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, @@ -73,65 +106,49 @@ class Exact_Alphacomplex_dD : public Abstract_alpha_complex { return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); } + virtual std::size_t num_vertices() const { + return alpha_complex_.num_vertices(); + } + private: bool exact_version_; - Alpha_complex<Kernel> alpha_complex_; + Alpha_complex<Kernel, Weighted> alpha_complex_; }; -class Inexact_Alphacomplex_dD : public Abstract_alpha_complex { +template <bool Weighted = false> +class Inexact_alpha_complex_dD final : public Abstract_alpha_complex { private: using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; - using Point = typename Kernel::Point_d; + using Bare_point = typename Kernel::Point_d; + using Point = std::conditional_t<Weighted, typename Kernel::Weighted_point_d, + typename Kernel::Point_d>; public: - Inexact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version) - : exact_version_(exact_version), - alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>)) { + Inexact_alpha_complex_dD(const std::vector<std::vector<double>>& points) + : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>)) { + } + + Inexact_alpha_complex_dD(const std::vector<std::vector<double>>& points, const std::vector<double>& weights) + : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>), weights) { } virtual std::vector<double> get_point(int vh) override { - Point const& point = alpha_complex_.get_point(vh); - return pt_cgal_to_cython(point); + // Can be a Weighted or a Bare point in function of Weighted + return Point_cgal_to_cython<Point, Weighted>()(alpha_complex_.get_point(vh)); } virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) override { - return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); + return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, false, default_filtration_value); } - private: - bool exact_version_; - Alpha_complex<Kernel> alpha_complex_; -}; - -template <complexity Complexity> -class Alphacomplex_3D : public Abstract_alpha_complex { - private: - using Point = typename Alpha_complex_3d<Complexity, false, false>::Bare_point_3; - - static Point pt_cython_to_cgal_3(std::vector<double> const& vec) { - return Point(vec[0], vec[1], vec[2]); - } - - public: - Alphacomplex_3D(const std::vector<std::vector<double>>& points) - : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3)) { - } - - virtual std::vector<double> get_point(int vh) override { - Point const& point = alpha_complex_.get_point(vh); - return pt_cgal_to_cython(point); - } - - virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, - bool default_filtration_value) override { - return alpha_complex_.create_complex(*simplex_tree, max_alpha_square); + virtual std::size_t num_vertices() const { + return alpha_complex_.num_vertices(); } private: - Alpha_complex_3d<Complexity, false, false> alpha_complex_; + Alpha_complex<Kernel, Weighted> alpha_complex_; }; - } // namespace alpha_complex } // namespace Gudhi diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 23be194d..671af4a4 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -27,10 +27,23 @@ namespace alpha_complex { class Alpha_complex_interface { public: - Alpha_complex_interface(const std::vector<std::vector<double>>& points, bool fast_version, bool exact_version) - : points_(points), - fast_version_(fast_version), - exact_version_(exact_version) { + Alpha_complex_interface(const std::vector<std::vector<double>>& points, + const std::vector<double>& weights, + bool fast_version, bool exact_version) { + const bool weighted = (weights.size() > 0); + if (fast_version) { + if (weighted) { + alpha_ptr_ = std::make_unique<Inexact_alpha_complex_dD<true>>(points, weights); + } else { + alpha_ptr_ = std::make_unique<Inexact_alpha_complex_dD<false>>(points); + } + } else { + if (weighted) { + alpha_ptr_ = std::make_unique<Exact_alpha_complex_dD<true>>(points, weights, exact_version); + } else { + alpha_ptr_ = std::make_unique<Exact_alpha_complex_dD<false>>(points, exact_version); + } + } } std::vector<double> get_point(int vh) { @@ -39,38 +52,13 @@ class Alpha_complex_interface { void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) { - if (points_.size() > 0) { - std::size_t dimension = points_[0].size(); - if (dimension == 3 && !default_filtration_value) { - if (fast_version_) - alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::FAST>>(points_); - else if (exact_version_) - alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::EXACT>>(points_); - else - alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::SAFE>>(points_); - if (!alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value)) { - // create_simplex_tree will fail if all points are on a plane - Retry with dD by setting dimension to 2 - dimension--; - alpha_ptr_.reset(); - } - } - // Not ** else ** because we have to take into account if 3d fails - if (dimension != 3 || default_filtration_value) { - if (fast_version_) { - alpha_ptr_ = std::make_unique<Inexact_Alphacomplex_dD>(points_, exact_version_); - } else { - alpha_ptr_ = std::make_unique<Exact_Alphacomplex_dD>(points_, exact_version_); - } - alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); - } - } + // Nothing to be done in case of an empty point set + if (alpha_ptr_->num_vertices() > 0) + alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); } private: std::unique_ptr<Abstract_alpha_complex> alpha_ptr_; - std::vector<std::vector<double>> points_; - bool fast_version_; - bool exact_version_; }; } // namespace alpha_complex diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 3061884f..95c0c037 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -15,10 +15,15 @@ #include <gudhi/distance_functions.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Points_off_io.h> +#ifdef GUDHI_USE_EIGEN3 +#include <gudhi/Flag_complex_edge_collapser.h> +#endif #include <iostream> #include <vector> #include <utility> // std::pair +#include <tuple> +#include <iterator> // for std::distance namespace Gudhi { @@ -36,6 +41,7 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator; using Complex_simplex_iterator = typename Base::Complex_simplex_iterator; using Extended_filtration_data = typename Base::Extended_filtration_data; + using Boundary_simplex_iterator = typename Base::Boundary_simplex_iterator; using Siblings = typename Base::Siblings; using Node = typename Base::Node; @@ -183,6 +189,38 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { return new_dgm; } + Simplex_tree_interface* collapse_edges(int nb_collapse_iteration) { +#ifdef GUDHI_USE_EIGEN3 + using Filtered_edge = std::tuple<Vertex_handle, Vertex_handle, Filtration_value>; + std::vector<Filtered_edge> edges; + for (Simplex_handle sh : Base::skeleton_simplex_range(1)) { + if (Base::dimension(sh) == 1) { + typename Base::Simplex_vertex_range rg = Base::simplex_vertex_range(sh); + auto vit = rg.begin(); + Vertex_handle v = *vit; + Vertex_handle w = *++vit; + edges.emplace_back(v, w, Base::filtration(sh)); + } + } + + for (int iteration = 0; iteration < nb_collapse_iteration; iteration++) { + edges = Gudhi::collapse::flag_complex_collapse_edges(edges); + } + Simplex_tree_interface* collapsed_stree_ptr = new Simplex_tree_interface(); + // Copy the original 0-skeleton + for (Simplex_handle sh : Base::skeleton_simplex_range(0)) { + collapsed_stree_ptr->insert({*(Base::simplex_vertex_range(sh).begin())}, Base::filtration(sh)); + } + // Insert remaining edges + for (auto remaining_edge : edges) { + collapsed_stree_ptr->insert({std::get<0>(remaining_edge), std::get<1>(remaining_edge)}, std::get<2>(remaining_edge)); + } + return collapsed_stree_ptr; +#else + throw std::runtime_error("Unable to collapse edges as it requires Eigen3 >= 3.1.0."); +#endif + } + // Iterator over the simplex tree Complex_simplex_iterator get_simplices_iterator_begin() { // this specific case works because the range is just a pair of iterators - won't work if range was a vector @@ -214,6 +252,15 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { // this specific case works because the range is just a pair of iterators - won't work if range was a vector return Base::skeleton_simplex_range(dimension).end(); } + + std::pair<Boundary_simplex_iterator, Boundary_simplex_iterator> get_boundary_iterators(const Simplex& simplex) { + auto bd_sh = Base::find(simplex); + if (bd_sh == Base::null_simplex()) + throw std::runtime_error("simplex not found - cannot find boundaries"); + // this specific case works because the range is just a pair of iterators - won't work if range was a vector + auto boundary_srange = Base::boundary_simplex_range(bd_sh); + return std::make_pair(boundary_srange.begin(), boundary_srange.end()); + } }; } // namespace Gudhi diff --git a/src/python/include/Subsampling_interface.h b/src/python/include/Subsampling_interface.h index cdda851f..6aee7231 100644 --- a/src/python/include/Subsampling_interface.h +++ b/src/python/include/Subsampling_interface.h @@ -11,6 +11,7 @@ #ifndef INCLUDE_SUBSAMPLING_INTERFACE_H_ #define INCLUDE_SUBSAMPLING_INTERFACE_H_ +#include <gudhi/distance_functions.h> #include <gudhi/choose_n_farthest_points.h> #include <gudhi/pick_n_random_points.h> #include <gudhi/sparsify_point_set.h> @@ -27,14 +28,13 @@ namespace subsampling { using Subsampling_dynamic_kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; using Subsampling_point_d = Subsampling_dynamic_kernel::Point_d; -using Subsampling_ft = Subsampling_dynamic_kernel::FT; // ------ choose_n_farthest_points ------ std::vector<std::vector<double>> subsampling_n_farthest_points(const std::vector<std::vector<double>>& points, unsigned nb_points) { std::vector<std::vector<double>> landmarks; - Subsampling_dynamic_kernel k; - choose_n_farthest_points(k, points, nb_points, random_starting_point, std::back_inserter(landmarks)); + choose_n_farthest_points(Euclidean_distance(), points, nb_points, + random_starting_point, std::back_inserter(landmarks)); return landmarks; } @@ -42,8 +42,8 @@ std::vector<std::vector<double>> subsampling_n_farthest_points(const std::vector std::vector<std::vector<double>> subsampling_n_farthest_points(const std::vector<std::vector<double>>& points, unsigned nb_points, unsigned starting_point) { std::vector<std::vector<double>> landmarks; - Subsampling_dynamic_kernel k; - choose_n_farthest_points(k, points, nb_points, starting_point, std::back_inserter(landmarks)); + choose_n_farthest_points(Euclidean_distance(), points, nb_points, + starting_point, std::back_inserter(landmarks)); return landmarks; } diff --git a/src/python/pyproject.toml b/src/python/pyproject.toml new file mode 100644 index 00000000..55b64466 --- /dev/null +++ b/src/python/pyproject.toml @@ -0,0 +1,3 @@ +[build-system] +requires = ["setuptools>=24.2.0", "wheel", "numpy>=1.15.0", "cython>=0.27", "pybind11"] +build-backend = "setuptools.build_meta" diff --git a/src/python/setup.py.in b/src/python/setup.py.in index 98d058fc..2c67c2c5 100644 --- a/src/python/setup.py.in +++ b/src/python/setup.py.in @@ -5,6 +5,7 @@ Copyright (C) 2019 Inria Modification(s): + - 2021/12 Vincent Rouvreau: Python 3.5 as minimal version - YYYY/MM Author: Description of the modification """ @@ -41,10 +42,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': '3'}) for module in pybind11_modules: my_include_dirs = include_dirs + [pybind11.get_include(False), pybind11.get_include(True)] @@ -72,7 +72,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 +83,11 @@ 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',], + python_requires='>=3.5.0', + install_requires = ['numpy >= 1.15.0',], package_data={"": ["*.dll"], }, ) diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py index a4ee260b..f15284f3 100755 --- a/src/python/test/test_alpha_complex.py +++ b/src/python/test/test_alpha_complex.py @@ -8,10 +8,12 @@ - YYYY/MM Author: Description of the modification """ -import gudhi as gd +from gudhi import AlphaComplex import math import numpy as np import pytest +import warnings + try: # python3 from itertools import zip_longest @@ -19,22 +21,24 @@ except ImportError: # python2 from itertools import izip_longest as zip_longest -__author__ = "Vincent Rouvreau" -__copyright__ = "Copyright (C) 2016 Inria" -__license__ = "MIT" def _empty_alpha(precision): - alpha_complex = gd.AlphaComplex(points=[[0, 0]], precision = precision) + alpha_complex = AlphaComplex(precision = precision) + assert alpha_complex.__is_defined() == True + +def _one_2d_point_alpha(precision): + alpha_complex = AlphaComplex(points=[[0, 0]], precision = precision) assert alpha_complex.__is_defined() == True def test_empty_alpha(): for precision in ['fast', 'safe', 'exact']: _empty_alpha(precision) + _one_2d_point_alpha(precision) def _infinite_alpha(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - alpha_complex = gd.AlphaComplex(points=point_list, precision = precision) + alpha_complex = AlphaComplex(points=point_list, precision = precision) assert alpha_complex.__is_defined() == True simplex_tree = alpha_complex.create_simplex_tree() @@ -69,18 +73,9 @@ def _infinite_alpha(precision): assert point_list[1] == alpha_complex.get_point(1) assert point_list[2] == alpha_complex.get_point(2) assert point_list[3] == alpha_complex.get_point(3) - try: - alpha_complex.get_point(4) == [] - except IndexError: - pass - else: - assert False - try: - alpha_complex.get_point(125) == [] - except IndexError: - pass - else: - assert False + + with pytest.raises(IndexError): + alpha_complex.get_point(len(point_list)) def test_infinite_alpha(): for precision in ['fast', 'safe', 'exact']: @@ -88,7 +83,7 @@ def test_infinite_alpha(): def _filtered_alpha(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - filtered_alpha = gd.AlphaComplex(points=point_list, precision = precision) + filtered_alpha = AlphaComplex(points=point_list, precision = precision) simplex_tree = filtered_alpha.create_simplex_tree(max_alpha_square=0.25) @@ -99,18 +94,9 @@ def _filtered_alpha(precision): assert point_list[1] == filtered_alpha.get_point(1) assert point_list[2] == filtered_alpha.get_point(2) assert point_list[3] == filtered_alpha.get_point(3) - try: - filtered_alpha.get_point(4) == [] - except IndexError: - pass - else: - assert False - try: - filtered_alpha.get_point(125) == [] - except IndexError: - pass - else: - assert False + + with pytest.raises(IndexError): + filtered_alpha.get_point(len(point_list)) assert list(simplex_tree.get_filtration()) == [ ([0], 0.0), @@ -141,10 +127,10 @@ def _safe_alpha_persistence_comparison(precision): embedding2 = [[signal[i], delayed[i]] for i in range(len(time))] #build alpha complex and simplex tree - alpha_complex1 = gd.AlphaComplex(points=embedding1, precision = precision) + alpha_complex1 = AlphaComplex(points=embedding1, precision = precision) simplex_tree1 = alpha_complex1.create_simplex_tree() - alpha_complex2 = gd.AlphaComplex(points=embedding2, precision = precision) + alpha_complex2 = AlphaComplex(points=embedding2, precision = precision) simplex_tree2 = alpha_complex2.create_simplex_tree() diag1 = simplex_tree1.persistence() @@ -162,7 +148,7 @@ def test_safe_alpha_persistence_comparison(): def _delaunay_complex(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - filtered_alpha = gd.AlphaComplex(points=point_list, precision = precision) + filtered_alpha = AlphaComplex(points=point_list, precision = precision) simplex_tree = filtered_alpha.create_simplex_tree(default_filtration_value = True) @@ -173,18 +159,11 @@ def _delaunay_complex(precision): assert point_list[1] == filtered_alpha.get_point(1) assert point_list[2] == filtered_alpha.get_point(2) assert point_list[3] == filtered_alpha.get_point(3) - try: - filtered_alpha.get_point(4) == [] - except IndexError: - pass - else: - assert False - try: - filtered_alpha.get_point(125) == [] - except IndexError: - pass - else: - assert False + + with pytest.raises(IndexError): + filtered_alpha.get_point(4) + with pytest.raises(IndexError): + filtered_alpha.get_point(125) for filtered_value in simplex_tree.get_filtration(): assert math.isnan(filtered_value[1]) @@ -198,8 +177,13 @@ def test_delaunay_complex(): _delaunay_complex(precision) def _3d_points_on_a_plane(precision, default_filtration_value): - alpha = gd.AlphaComplex(off_file=gd.__root_source_dir__ + '/data/points/alphacomplexdoc.off', - precision = precision) + alpha = AlphaComplex(points = [[1.0, 1.0 , 0.0], + [7.0, 0.0 , 0.0], + [4.0, 6.0 , 0.0], + [9.0, 6.0 , 0.0], + [0.0, 14.0, 0.0], + [2.0, 19.0, 0.0], + [9.0, 17.0, 0.0]], precision = precision) simplex_tree = alpha.create_simplex_tree(default_filtration_value = default_filtration_value) assert simplex_tree.dimension() == 2 @@ -213,10 +197,10 @@ def test_3d_points_on_a_plane(): def _3d_tetrahedrons(precision): points = 10*np.random.rand(10, 3) - alpha = gd.AlphaComplex(points=points, precision = precision) + alpha = AlphaComplex(points = points, precision = precision) st_alpha = alpha.create_simplex_tree(default_filtration_value = False) # New AlphaComplex for get_point to work - delaunay = gd.AlphaComplex(points=points, precision = precision) + delaunay = AlphaComplex(points = points, precision = precision) st_delaunay = delaunay.create_simplex_tree(default_filtration_value = True) delaunay_tetra = [] @@ -245,3 +229,60 @@ def _3d_tetrahedrons(precision): def test_3d_tetrahedrons(): for precision in ['fast', 'safe', 'exact']: _3d_tetrahedrons(precision) + +def test_off_file_deprecation_warning(): + off_file = open("alphacomplexdoc.off", "w") + off_file.write("OFF \n" \ + "7 0 0 \n" \ + "1.0 1.0 0.0\n" \ + "7.0 0.0 0.0\n" \ + "4.0 6.0 0.0\n" \ + "9.0 6.0 0.0\n" \ + "0.0 14.0 0.0\n" \ + "2.0 19.0 0.0\n" \ + "9.0 17.0 0.0\n" ) + off_file.close() + + with pytest.warns(DeprecationWarning): + alpha = AlphaComplex(off_file="alphacomplexdoc.off") + +def test_non_existing_off_file(): + with pytest.warns(DeprecationWarning): + with pytest.raises(FileNotFoundError): + alpha = AlphaComplex(off_file="pouetpouettralala.toubiloubabdou") + +def test_inconsistency_points_and_weights(): + points = [[1.0, 1.0 , 0.0], + [7.0, 0.0 , 0.0], + [4.0, 6.0 , 0.0], + [9.0, 6.0 , 0.0], + [0.0, 14.0, 0.0], + [2.0, 19.0, 0.0], + [9.0, 17.0, 0.0]] + with pytest.raises(ValueError): + # 7 points, 8 weights, on purpose + alpha = AlphaComplex(points = points, + weights = [1., 2., 3., 4., 5., 6., 7., 8.]) + + with pytest.raises(ValueError): + # 7 points, 6 weights, on purpose + alpha = AlphaComplex(points = points, + weights = [1., 2., 3., 4., 5., 6.]) + +def _weighted_doc_example(precision): + stree = AlphaComplex(points=[[ 1., -1., -1.], + [-1., 1., -1.], + [-1., -1., 1.], + [ 1., 1., 1.], + [ 2., 2., 2.]], + weights = [4., 4., 4., 4., 1.], + precision = precision).create_simplex_tree() + + assert stree.filtration([0, 1, 2, 3]) == pytest.approx(-1.) + assert stree.filtration([0, 1, 3, 4]) == pytest.approx(95.) + assert stree.filtration([0, 2, 3, 4]) == pytest.approx(95.) + assert stree.filtration([1, 2, 3, 4]) == pytest.approx(95.) + +def test_weighted_doc_example(): + for precision in ['fast', 'safe', 'exact']: + _weighted_doc_example(precision) diff --git a/src/python/test/test_betti_curve_representations.py b/src/python/test/test_betti_curve_representations.py new file mode 100755 index 00000000..6a45da4d --- /dev/null +++ b/src/python/test/test_betti_curve_representations.py @@ -0,0 +1,59 @@ +import numpy as np +import scipy.interpolate +import pytest + +from gudhi.representations.vector_methods import BettiCurve + +def test_betti_curve_is_irregular_betti_curve_followed_by_interpolation(): + m = 10 + n = 1000 + pinf = 0.05 + pzero = 0.05 + res = 100 + + pds = [] + for i in range(0, m): + pd = np.zeros((n, 2)) + pd[:, 0] = np.random.uniform(0, 10, n) + pd[:, 1] = np.random.uniform(pd[:, 0], 10, n) + pd[np.random.uniform(0, 1, n) < pzero, 0] = 0 + pd[np.random.uniform(0, 1, n) < pinf, 1] = np.inf + pds.append(pd) + + bc = BettiCurve(resolution=None, predefined_grid=None) + bc.fit(pds) + bettis = bc.transform(pds) + + bc2 = BettiCurve(resolution=None, predefined_grid=None) + bettis2 = bc2.fit_transform(pds) + assert((bc2.grid_ == bc.grid_).all()) + assert((bettis2 == bettis).all()) + + for i in range(0, m): + grid = np.linspace(pds[i][np.isfinite(pds[i])].min(), pds[i][np.isfinite(pds[i])].max() + 1, res) + bc_gridded = BettiCurve(predefined_grid=grid) + bc_gridded.fit([]) + bettis_gridded = bc_gridded(pds[i]) + + interp = scipy.interpolate.interp1d(bc.grid_, bettis[i, :], kind="previous", fill_value="extrapolate") + bettis_interp = np.array(interp(grid), dtype=int) + assert((bettis_interp == bettis_gridded).all()) + + +def test_empty_with_predefined_grid(): + random_grid = np.sort(np.random.uniform(0, 1, 100)) + bc = BettiCurve(predefined_grid=random_grid) + bettis = bc.fit_transform([]) + assert((bc.grid_ == random_grid).all()) + assert((bettis == 0).all()) + + +def test_empty(): + bc = BettiCurve(resolution=None, predefined_grid=None) + bettis = bc.fit_transform([]) + assert(bc.grid_ == [-np.inf]) + assert((bettis == 0).all()) + +def test_wrong_value_of_predefined_grid(): + with pytest.raises(ValueError): + BettiCurve(predefined_grid=[1, 2, 3]) diff --git a/src/python/test/test_bottleneck_distance.py b/src/python/test/test_bottleneck_distance.py index 6915bea8..07fcc9cc 100755 --- a/src/python/test/test_bottleneck_distance.py +++ b/src/python/test/test_bottleneck_distance.py @@ -25,3 +25,15 @@ def test_basic_bottleneck(): assert gudhi.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, abs=0.1) assert gudhi.hera.bottleneck_distance(diag1, diag2, 0) == 0.75 assert gudhi.hera.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, rel=0.1) + + import numpy as np + + # Translating both diagrams along the diagonal should not affect the distance, + # checks that negative numbers are not an issue + diag1 = np.array(diag1) - 100 + diag2 = np.array(diag2) - 100 + + assert gudhi.bottleneck_distance(diag1, diag2) == 0.75 + assert gudhi.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, abs=0.1) + assert gudhi.hera.bottleneck_distance(diag1, diag2, 0) == 0.75 + assert gudhi.hera.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, rel=0.1) 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..fdfddc4b 100755 --- a/src/python/test/test_reader_utils.py +++ b/src/python/test/test_reader_utils.py @@ -8,8 +8,9 @@ - YYYY/MM Author: Description of the modification """ -import gudhi +import gudhi as gd import numpy as np +from pytest import raises __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2017 Inria" @@ -18,7 +19,7 @@ __license__ = "MIT" def test_non_existing_csv_file(): # Try to open a non existing file - matrix = gudhi.read_lower_triangular_matrix_from_csv_file( + matrix = gd.read_lower_triangular_matrix_from_csv_file( csv_file="pouetpouettralala.toubiloubabdou" ) assert matrix == [] @@ -29,8 +30,8 @@ def test_full_square_distance_matrix_csv_file(): test_file = open("full_square_distance_matrix.csv", "w") 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" + matrix = gd.read_lower_triangular_matrix_from_csv_file( + csv_file="full_square_distance_matrix.csv", separator=";" ) assert matrix == [[], [1.0], [2.0, 4.0], [3.0, 5.0, 6.0]] @@ -40,7 +41,7 @@ def test_lower_triangular_distance_matrix_csv_file(): test_file = open("lower_triangular_distance_matrix.csv", "w") test_file.write("\n1,\n2,3,\n4,5,6,\n7,8,9,10,") test_file.close() - matrix = gudhi.read_lower_triangular_matrix_from_csv_file( + matrix = gd.read_lower_triangular_matrix_from_csv_file( csv_file="lower_triangular_distance_matrix.csv", separator="," ) assert matrix == [[], [1.0], [2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0, 10.0]] @@ -48,11 +49,11 @@ def test_lower_triangular_distance_matrix_csv_file(): def test_non_existing_persistence_file(): # Try to open a non existing file - persistence = gudhi.read_persistence_intervals_grouped_by_dimension( + persistence = gd.read_persistence_intervals_grouped_by_dimension( persistence_file="pouetpouettralala.toubiloubabdou" ) assert persistence == [] - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="pouetpouettralala.toubiloubabdou", only_this_dim=1 ) np.testing.assert_array_equal(persistence, []) @@ -65,21 +66,21 @@ def test_read_persistence_intervals_without_dimension(): "# Simple persistence diagram without dimension\n2.7 3.7\n9.6 14.\n34.2 34.974\n3. inf" ) test_file.close() - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_without_dimension.pers" ) np.testing.assert_array_equal( persistence, [(2.7, 3.7), (9.6, 14.0), (34.2, 34.974), (3.0, float("Inf"))] ) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_without_dimension.pers", only_this_dim=0 ) np.testing.assert_array_equal(persistence, []) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_without_dimension.pers", only_this_dim=1 ) np.testing.assert_array_equal(persistence, []) - persistence = gudhi.read_persistence_intervals_grouped_by_dimension( + persistence = gd.read_persistence_intervals_grouped_by_dimension( persistence_file="persistence_intervals_without_dimension.pers" ) assert persistence == { @@ -94,29 +95,29 @@ def test_read_persistence_intervals_with_dimension(): "# Simple persistence diagram with dimension\n0 2.7 3.7\n1 9.6 14.\n3 34.2 34.974\n1 3. inf" ) test_file.close() - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_with_dimension.pers" ) np.testing.assert_array_equal( persistence, [(2.7, 3.7), (9.6, 14.0), (34.2, 34.974), (3.0, float("Inf"))] ) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=0 ) np.testing.assert_array_equal(persistence, [(2.7, 3.7)]) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=1 ) np.testing.assert_array_equal(persistence, [(9.6, 14.0), (3.0, float("Inf"))]) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=2 ) np.testing.assert_array_equal(persistence, []) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=3 ) np.testing.assert_array_equal(persistence, [(34.2, 34.974)]) - persistence = gudhi.read_persistence_intervals_grouped_by_dimension( + persistence = gd.read_persistence_intervals_grouped_by_dimension( persistence_file="persistence_intervals_with_dimension.pers" ) assert persistence == { diff --git a/src/python/test/test_representations.py b/src/python/test/test_representations.py index 589cee00..d219ce7a 100755 --- a/src/python/test/test_representations.py +++ b/src/python/test/test_representations.py @@ -3,6 +3,22 @@ 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(): @@ -15,6 +31,7 @@ def test_representations_examples(): return None +from gudhi.representations.vector_methods import Atol from gudhi.representations.metrics import * from gudhi.representations.kernel_methods import * @@ -36,8 +53,118 @@ def test_multiple(): d2 = BottleneckDistance(epsilon=0.00001).fit_transform(l1) d3 = pairwise_persistence_diagram_distances(l1, l1b, e=0.00001, n_jobs=4) assert d1 == pytest.approx(d2) - assert d3 == pytest.approx(d2, abs=1e-5) # Because of 0 entries (on the diagonal) + assert d3 == pytest.approx(d2, abs=1e-5) # Because of 0 entries (on the diagonal) d1 = pairwise_persistence_diagram_distances(l1, l2, metric="wasserstein", order=2, internal_p=2) d2 = WassersteinDistance(order=2, internal_p=2, n_jobs=4).fit(l2).transform(l1) print(d1.shape, d2.shape) - assert d1 == pytest.approx(d2, rel=.02) + assert d1 == pytest.approx(d2, rel=0.02) + + +# Test sorted values as points order can be inverted, and sorted test is not documentation-friendly +# Note the test below must be up to date with the Atol class documentation +def test_atol_doc(): + a = np.array([[1, 2, 4], [1, 4, 0], [1, 0, 4]]) + b = np.array([[4, 2, 0], [4, 4, 0], [4, 0, 2]]) + c = np.array([[3, 2, -1], [1, 2, -1]]) + + atol_vectoriser = Atol(quantiser=KMeans(n_clusters=2, random_state=202006)) + # Atol will do + # X = np.concatenate([a,b,c]) + # kmeans = KMeans(n_clusters=2, random_state=202006).fit(X) + # kmeans.labels_ will be : array([1, 0, 1, 0, 0, 1, 0, 0]) + first_cluster = np.asarray([a[0], a[2], b[2]]) + second_cluster = np.asarray([a[1], b[0], b[2], c[0], c[1]]) + + # Check the center of the first_cluster and second_cluster are in Atol centers + centers = atol_vectoriser.fit(X=[a, b, c]).centers + np.isclose(centers, first_cluster.mean(axis=0)).all(1).any() + np.isclose(centers, second_cluster.mean(axis=0)).all(1).any() + + vectorization = atol_vectoriser.transform(X=[a, b, c]) + assert np.allclose(vectorization[0], atol_vectoriser(a)) + assert np.allclose(vectorization[1], atol_vectoriser(b)) + assert np.allclose(vectorization[2], atol_vectoriser(c)) + + +def test_dummy_atol(): + a = np.array([[1, 2, 4], [1, 4, 0], [1, 0, 4]]) + b = np.array([[4, 2, 0], [4, 4, 0], [4, 0, 2]]) + c = np.array([[3, 2, -1], [1, 2, -1]]) + + for weighting_method in ["cloud", "iidproba"]: + for contrast in ["gaussian", "laplacian", "indicator"]: + atol_vectoriser = Atol( + quantiser=KMeans(n_clusters=1, random_state=202006), + weighting_method=weighting_method, + contrast=contrast, + ) + atol_vectoriser.fit([a, b, c]) + atol_vectoriser(a) + atol_vectoriser.transform(X=[a, b, c]) + + +from gudhi.representations.vector_methods import BettiCurve + +def test_infinity(): + a = np.array([[1.0, 8.0], [2.0, np.inf], [3.0, 4.0]]) + c = BettiCurve(20, [0.0, 10.0])(a) + 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 34173e78..0436d891 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -8,7 +8,8 @@ - YYYY/MM Author: Description of the modification """ -from gudhi import SimplexTree +from gudhi import SimplexTree, __GUDHI_USE_EIGEN3 +import numpy as np import pytest import numpy as np @@ -344,13 +345,185 @@ def test_simplices_iterator(): print("filtration is: ", simplex[1]) assert st.filtration(simplex[0]) == simplex[1] +def test_collapse_edges(): + st = SimplexTree() + + assert st.insert([0, 1], filtration=1.0) == True + assert st.insert([1, 2], filtration=1.0) == True + assert st.insert([2, 3], filtration=1.0) == True + assert st.insert([0, 3], filtration=1.0) == True + assert st.insert([0, 2], filtration=2.0) == True + assert st.insert([1, 3], filtration=2.0) == True + + assert st.num_simplices() == 10 + + if __GUDHI_USE_EIGEN3: + st.collapse_edges() + assert st.num_simplices() == 9 + assert st.find([1, 3]) == False + for simplex in st.get_skeleton(0): + assert simplex[1] == 1. + else: + # If no Eigen3, collapse_edges throws an exception + with pytest.raises(RuntimeError): + st.collapse_edges() + +def test_reset_filtration(): + st = SimplexTree() + + assert st.insert([0, 1, 2], 3.) == True + assert st.insert([0, 3], 2.) == True + assert st.insert([3, 4, 5], 3.) == True + assert st.insert([0, 1, 6, 7], 4.) == True + + # Guaranteed by construction + for simplex in st.get_simplices(): + assert st.filtration(simplex[0]) >= 2. + + # dimension until 5 even if simplex tree is of dimension 3 to test the limits + for dimension in range(5, -1, -1): + st.reset_filtration(0., dimension) + for simplex in st.get_skeleton(3): + print(simplex) + if len(simplex[0]) < (dimension) + 1: + assert st.filtration(simplex[0]) >= 2. + else: + assert st.filtration(simplex[0]) == 0. + +def test_boundaries_iterator(): + st = SimplexTree() + + assert st.insert([0, 1, 2, 3], filtration=1.0) == True + assert st.insert([1, 2, 3, 4], filtration=2.0) == True + + assert list(st.get_boundaries([1, 2, 3])) == [([1, 2], 1.0), ([1, 3], 1.0), ([2, 3], 1.0)] + assert list(st.get_boundaries([2, 3, 4])) == [([2, 3], 1.0), ([2, 4], 2.0), ([3, 4], 2.0)] + assert list(st.get_boundaries([2])) == [] + + with pytest.raises(RuntimeError): + list(st.get_boundaries([])) + + with pytest.raises(RuntimeError): + list(st.get_boundaries([0, 4])) # (0, 4) does not exist + + 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() + + assert st1 == st2 + + st1.insert([1,2,3], 4.) + assert st1 != st2 + + st2.insert([1,2,3], 4.) + assert st1 == st2 + +def test_simplex_tree_deep_copy(): + st = SimplexTree() + st.insert([1, 2, 3], 0.) + # compute persistence only on the original + st.compute_persistence() + + st_copy = st.copy() + assert st_copy == st + st_filt_list = list(st.get_filtration()) + + # check persistence is not copied + assert st.__is_persistence_defined() == True + assert st_copy.__is_persistence_defined() == False + + # remove something in the copy and check the copy is included in the original + st_copy.remove_maximal_simplex([1, 2, 3]) + a_filt_list = list(st_copy.get_filtration()) + assert len(a_filt_list) < len(st_filt_list) + + for a_splx in a_filt_list: + assert a_splx in st_filt_list + + # test double free + del st + del st_copy + +def test_simplex_tree_deep_copy_constructor(): + st = SimplexTree() + st.insert([1, 2, 3], 0.) + # compute persistence only on the original + st.compute_persistence() + + st_copy = SimplexTree(st) + assert st_copy == st + st_filt_list = list(st.get_filtration()) + + # check persistence is not copied + assert st.__is_persistence_defined() == True + assert st_copy.__is_persistence_defined() == False + + # remove something in the copy and check the copy is included in the original + st_copy.remove_maximal_simplex([1, 2, 3]) + a_filt_list = list(st_copy.get_filtration()) + assert len(a_filt_list) < len(st_filt_list) + + for a_splx in a_filt_list: + assert a_splx in st_filt_list + + # test double free + del st + del st_copy + +def test_simplex_tree_constructor_exception(): + with pytest.raises(TypeError): + st = SimplexTree(other = "Construction from a string shall raise an exception") def test_create_from_array(): a = np.array([[1, 4, 13, 6], [4, 3, 11, 5], [13, 11, 10, 12], [6, 5, 12, 2]]) st = SimplexTree.create_from_array(a, max_filtration=5) assert list(st.get_filtration()) == [([0], 1.0), ([3], 2.0), ([1], 3.0), ([0, 1], 4.0), ([1, 3], 5.0)] - def test_insert_edges_from_coo_matrix(): try: from scipy.sparse import coo_matrix diff --git a/src/python/test/test_subsampling.py b/src/python/test/test_subsampling.py index 31f64e32..4019852e 100755 --- a/src/python/test/test_subsampling.py +++ b/src/python/test/test_subsampling.py @@ -141,12 +141,16 @@ def test_simple_sparsify_points(): # assert gudhi.sparsify_point_set(points = [], min_squared_dist = 0.0) == [] # assert gudhi.sparsify_point_set(points = [], min_squared_dist = 10.0) == [] assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=0.0) == point_set - assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=1.0) == point_set - assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=2.0) == [ + assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=0.999) == point_set + assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=1.001) == [ [0, 1], [1, 0], ] - assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=2.01) == [[0, 1]] + assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=1.999) == [ + [0, 1], + [1, 0], + ] + assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=2.001) == [[0, 1]] assert ( len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=0.0)) @@ -157,11 +161,11 @@ def test_simple_sparsify_points(): == 5 ) assert ( - len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=40.0)) + len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=40.1)) == 4 ) assert ( - len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=90.0)) + len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=89.9)) == 3 ) assert ( @@ -169,7 +173,7 @@ def test_simple_sparsify_points(): == 2 ) assert ( - len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=325.0)) + len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=324.9)) == 2 ) assert ( 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 90d26809..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,39 +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) -def test_wasserstein_distance_grad(): - import torch - - diag1 = torch.tensor([[2.7, 3.7], [9.6, 14.0], [34.2, 34.974]], requires_grad=True) - diag2 = torch.tensor([[2.8, 4.45], [9.5, 14.1]], requires_grad=True) - diag3 = torch.tensor([[2.8, 4.45], [9.5, 14.1]], requires_grad=True) - assert diag1.grad is None and diag2.grad is None and diag3.grad is None - dist12 = pot(diag1, diag2, internal_p=2, order=2, enable_autodiff=True) - dist30 = pot(diag3, torch.tensor([]), internal_p=2, order=2, enable_autodiff=True) - dist12.backward() - dist30.backward() - assert not torch.isnan(diag1.grad).any() and not torch.isnan(diag2.grad).any() and not torch.isnan(diag3.grad).any() - diag4 = torch.tensor([[0., 10.]], requires_grad=True) - diag5 = torch.tensor([[1., 11.], [3., 4.]], requires_grad=True) - dist45 = pot(diag4, diag5, internal_p=1, order=1, enable_autodiff=True) - assert dist45 == 3. - dist45.backward() - assert np.array_equal(diag4.grad, [[-1., -1.]]) - assert np.array_equal(diag5.grad, [[1., 1.], [-1., 1.]]) - diag6 = torch.tensor([[5., 10.]], requires_grad=True) - pot(diag6, diag6, internal_p=2, order=2, enable_autodiff=True).backward() - # https://github.com/jonasrauber/eagerpy/issues/6 - # assert np.array_equal(diag6.grad, [[0., 0.]]) diff --git a/src/python/test/test_wasserstein_with_tensors.py b/src/python/test/test_wasserstein_with_tensors.py new file mode 100755 index 00000000..e3f1411a --- /dev/null +++ b/src/python/test/test_wasserstein_with_tensors.py @@ -0,0 +1,47 @@ +""" 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): Mathieu Carriere + + Copyright (C) 2020 Inria + + Modification(s): + - YYYY/MM Author: Description of the modification +""" + +from gudhi.wasserstein import wasserstein_distance as pot +import numpy as np +import torch +import tensorflow as tf + +def test_wasserstein_distance_grad(): + diag1 = torch.tensor([[2.7, 3.7], [9.6, 14.0], [34.2, 34.974]], requires_grad=True) + diag2 = torch.tensor([[2.8, 4.45], [9.5, 14.1]], requires_grad=True) + diag3 = torch.tensor([[2.8, 4.45], [9.5, 14.1]], requires_grad=True) + assert diag1.grad is None and diag2.grad is None and diag3.grad is None + dist12 = pot(diag1, diag2, internal_p=2, order=2, enable_autodiff=True) + dist30 = pot(diag3, torch.tensor([]), internal_p=2, order=2, enable_autodiff=True) + dist12.backward() + dist30.backward() + assert not torch.isnan(diag1.grad).any() and not torch.isnan(diag2.grad).any() and not torch.isnan(diag3.grad).any() + diag4 = torch.tensor([[0., 10.]], requires_grad=True) + diag5 = torch.tensor([[1., 11.], [3., 4.]], requires_grad=True) + dist45 = pot(diag4, diag5, internal_p=1, order=1, enable_autodiff=True) + assert dist45 == 3. + dist45.backward() + assert np.array_equal(diag4.grad, [[-1., -1.]]) + assert np.array_equal(diag5.grad, [[1., 1.], [-1., 1.]]) + diag6 = torch.tensor([[5., 10.]], requires_grad=True) + pot(diag6, diag6, internal_p=2, order=2, enable_autodiff=True).backward() + # https://github.com/jonasrauber/eagerpy/issues/6 + # assert np.array_equal(diag6.grad, [[0., 0.]]) + +def test_wasserstein_distance_grad_tensorflow(): + with tf.GradientTape() as tape: + diag4 = tf.convert_to_tensor(tf.Variable(initial_value=np.array([[0., 10.]]), trainable=True)) + diag5 = tf.convert_to_tensor(tf.Variable(initial_value=np.array([[1., 11.], [3., 4.]]), trainable=True)) + dist45 = pot(diag4, diag5, internal_p=1, order=1, enable_autodiff=True) + assert dist45 == 3. + + grads = tape.gradient(dist45, [diag4, diag5]) + assert np.array_equal(grads[0].values, [[-1., -1.]]) + assert np.array_equal(grads[1].values, [[1., 1.], [-1., 1.]])
\ No newline at end of file |
