diff options
Diffstat (limited to 'src/python')
27 files changed, 755 insertions, 532 deletions
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 5f323935..39e2acd4 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -44,7 +44,7 @@ function( add_gudhi_debug_info DEBUG_INFO ) endfunction( add_gudhi_debug_info ) if(PYTHONINTERP_FOUND) - if(PYBIND11_FOUND AND CYTHON_FOUND) + if(NUMPY_FOUND AND PYBIND11_FOUND AND CYTHON_FOUND) add_gudhi_debug_info("Pybind11 version ${PYBIND11_VERSION}") # PyBind11 modules set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'bottleneck', ") @@ -53,7 +53,7 @@ if(PYTHONINTERP_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}'off_utils', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'simplex_tree', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'rips_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'cubical_complex', ") @@ -129,9 +129,10 @@ if(PYTHONINTERP_FOUND) # Gudhi and CGAL compilation option if(MSVC) + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'/std:c++17', ") set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'/fp:strict', ") else(MSVC) - set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-std=c++14', ") + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-std=c++17', ") endif(MSVC) if(CMAKE_COMPILER_IS_GNUCXX) set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-frounding-math', ") @@ -151,7 +152,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_EIGEN3_ENABLED', ") endif (EIGEN3_FOUND) - set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_reader', ") + set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_utils', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'simplex_tree', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'rips_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'cubical_complex', ") @@ -162,10 +163,10 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'clustering/_tomato', ") set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera/wasserstein', ") set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera/bottleneck', ") + set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'nerve_gic', ") 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}'subsampling', ") @@ -246,8 +247,8 @@ if(PYTHONINTERP_FOUND) # Specific for Mac if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin") - set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-mmacosx-version-min=10.12', ") - set(GUDHI_PYTHON_EXTRA_LINK_ARGS "${GUDHI_PYTHON_EXTRA_LINK_ARGS}'-mmacosx-version-min=10.12', ") + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-mmacosx-version-min=10.14', ") + set(GUDHI_PYTHON_EXTRA_LINK_ARGS "${GUDHI_PYTHON_EXTRA_LINK_ARGS}'-mmacosx-version-min=10.14', ") endif(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") # Loop on INCLUDE_DIRECTORIES PROPERTY @@ -431,38 +432,38 @@ if(PYTHONINTERP_FOUND) ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/bottleneck_basic_example.py") add_gudhi_py_test(test_bottleneck_distance) + endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) - # Cover complex - file(COPY ${CMAKE_SOURCE_DIR}/data/points/human.off DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) - file(COPY ${CMAKE_SOURCE_DIR}/data/points/COIL_database/lucky_cat.off DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) - 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 "${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) + # Cover complex + file(COPY ${CMAKE_SOURCE_DIR}/data/points/human.off DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + file(COPY ${CMAKE_SOURCE_DIR}/data/points/COIL_database/lucky_cat.off DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + 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 "${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 "${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_coordinate_gic_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/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 "${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_functional_gic_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/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 "${GUDHI_PYTHON_PATH_ENV}" - ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/voronoi_graph_induced_complex.py" - -f human.off -n 700 -v) + add_test(NAME cover_complex_voronoi_gic_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/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) + add_gudhi_py_test(test_cover_complex) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) # Alpha @@ -545,6 +546,7 @@ if(PYTHONINTERP_FOUND) # Reader utils add_gudhi_py_test(test_reader_utils) + add_gudhi_py_test(test_off) # Wasserstein if(OT_FOUND) @@ -621,10 +623,10 @@ if(PYTHONINTERP_FOUND) # Set missing or not modules set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES") - else(PYBIND11_FOUND AND CYTHON_FOUND) - message("++ Python module will not be compiled because cython and/or pybind11 was/were not found") + else(NUMPY_FOUND AND PYBIND11_FOUND AND CYTHON_FOUND) + message("++ Python module will not be compiled because numpy and/or cython and/or pybind11 was/were not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(PYBIND11_FOUND AND CYTHON_FOUND) + endif(NUMPY_FOUND AND PYBIND11_FOUND AND CYTHON_FOUND) else(PYTHONINTERP_FOUND) message("++ Python module will not be compiled because no Python interpreter was found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python" CACHE INTERNAL "GUDHI_MISSING_MODULES") diff --git a/src/python/doc/clustering.rst b/src/python/doc/clustering.rst index c5a57d3c..62422682 100644 --- a/src/python/doc/clustering.rst +++ b/src/python/doc/clustering.rst @@ -17,9 +17,8 @@ As a by-product, we produce the persistence diagram of the merge tree of the ini :include-source: import gudhi - f = open(gudhi.__root_source_dir__ + '/data/points/spiral_2d.csv', 'r') - import numpy as np - data = np.loadtxt(f) + from gudhi.datasets.remote import fetch_spiral_2d + data = fetch_spiral_2d() import matplotlib.pyplot as plt plt.scatter(data[:,0],data[:,1],marker='.',s=1) plt.show() diff --git a/src/python/doc/cubical_complex_sklearn_itf_ref.rst b/src/python/doc/cubical_complex_sklearn_itf_ref.rst index 2fb8ec6a..90ae9ccd 100644 --- a/src/python/doc/cubical_complex_sklearn_itf_ref.rst +++ b/src/python/doc/cubical_complex_sklearn_itf_ref.rst @@ -54,9 +54,9 @@ two holes in :math:`\mathbf{H}_1`, or, like in this example, three connected com X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=0) pipe = Pipeline( [ - ("cub_pers", CubicalPersistence(homology_dimensions=0, newshape=[28, 28], n_jobs=-2)), + ("cub_pers", CubicalPersistence(homology_dimensions=0, newshape=[-1, 28, 28], n_jobs=-2)), # Or for multiple persistence dimension computation - # ("cub_pers", CubicalPersistence(homology_dimensions=[0, 1], newshape=[28, 28], n_jobs=-2)), + # ("cub_pers", CubicalPersistence(homology_dimensions=[0, 1], newshape=[-1, 28, 28])), # ("H0_diags", DimensionSelector(index=0), # where index is the index in homology_dimensions array ("finite_diags", DiagramSelector(use=True, point_type="finite")), ( diff --git a/src/python/doc/installation.rst b/src/python/doc/installation.rst index 4eefd415..5491542f 100644 --- a/src/python/doc/installation.rst +++ b/src/python/doc/installation.rst @@ -39,7 +39,7 @@ 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.66.0, +The library uses c++17 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. @@ -150,7 +150,7 @@ You shall have something like: 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; + + Installed modules are: off_utils;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; @@ -188,7 +188,7 @@ A complete configuration would be : 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; + + Installed modules are: bottleneck;off_utils;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: @@ -391,7 +391,7 @@ The :doc:`persistence graphical tools </persistence_graphical_tools_user>` and 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. +`SciPy <http://scipy.org>`_ :math:`\geq` 1.6.0 as a backend if explicitly requested. TensorFlow ---------- diff --git a/src/python/doc/point_cloud.rst b/src/python/doc/point_cloud.rst index ffd8f85b..473b303f 100644 --- a/src/python/doc/point_cloud.rst +++ b/src/python/doc/point_cloud.rst @@ -13,6 +13,11 @@ File Readers .. autofunction:: gudhi.read_lower_triangular_matrix_from_csv_file +File Writers +------------ + +.. autofunction:: gudhi.write_points_to_off_file + Subsampling ----------- diff --git a/src/python/doc/representations_sum.inc b/src/python/doc/representations_sum.inc index 4298aea9..9515f044 100644 --- a/src/python/doc/representations_sum.inc +++ b/src/python/doc/representations_sum.inc @@ -1,14 +1,14 @@ .. table:: :widths: 30 40 30 - +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------+ - | .. 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 | - | | | | - | | | :License: MIT | - | | | | - | | | :Requires: `Scikit-learn <installation.html#scikit-learn>`_ | - +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------+ - | * :doc:`representations` | - +------------------------------------------------------------------+------------------------------------------------------------------------------------------------------------------------------+ + +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------------------+ + | .. figure:: | Vectorizations, distances and kernels that work on persistence | :Author: Mathieu Carrière, Martin Royer, Gard Spreemann, Wojciech Reise | + | img/sklearn-tda.png | diagrams, compatible with scikit-learn. | | + | | | :Since: GUDHI 3.1.0 | + | | | | + | | | :License: MIT | + | | | | + | | | :Requires: `Scikit-learn <installation.html#scikit-learn>`_ | + +------------------------------------------------------------------+----------------------------------------------------------------+-------------------------------------------------------------------------+ + | * :doc:`representations` | + +------------------------------------------------------------------+------------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/doc/rips_complex_user.rst b/src/python/doc/rips_complex_user.rst index c41a7803..a4e83462 100644 --- a/src/python/doc/rips_complex_user.rst +++ b/src/python/doc/rips_complex_user.rst @@ -34,9 +34,6 @@ A vertex name corresponds to the index of the point in the given range (aka. the On this example, as edges (4,5), (4,6) and (5,6) are in the complex, simplex (4,5,6) is added with the filtration value set with :math:`max(filtration(4,5), filtration(4,6), filtration(5,6))`. And so on for simplex (0,1,2,3). -If the :doc:`RipsComplex <rips_complex_ref>` interfaces are not detailed enough for your need, please refer to -rips_persistence_step_by_step.cpp C++ example, where the graph construction over the Simplex_tree is more detailed. - A Rips complex can easily become huge, even if we limit the length of the edges and the dimension of the simplices. One easy trick, before building a Rips complex on a point cloud, is to call :func:`~gudhi.sparsify_point_set` which removes points @@ -55,6 +52,13 @@ construction of a :class:`~gudhi.RipsComplex` object asks it to build a sparse R parameter :math:`\varepsilon=0.3`, while the default `sparse=None` builds the regular Rips complex. +Another option which is especially useful if you want to compute persistent homology in "high" dimension (2 or more, +sometimes even 1), is to build the Rips complex only up to dimension 1 (a graph), then use +:func:`~gudhi.SimplexTree.collapse_edges` to reduce the size of this graph, and finally call +:func:`~gudhi.SimplexTree.expansion` to get a simplicial complex of a suitable dimension to compute its homology. This +trick gives the same persistence diagram as one would get with a plain use of `RipsComplex`, with a complex that is +often significantly smaller and thus faster to process. + Point cloud ----------- @@ -117,54 +121,44 @@ Notice that if we use asking for a very sparse version (theory only gives some guarantee on the meaning of the output if `sparse<1`), 2 to 5 edges disappear, depending on the random vertex used to start the sparsification. -Example from OFF file -^^^^^^^^^^^^^^^^^^^^^ +Example step by step +^^^^^^^^^^^^^^^^^^^^ -This example builds the :doc:`RipsComplex <rips_complex_ref>` from the given -points in an OFF file, and max_edge_length value. -Then it creates a :doc:`SimplexTree <simplex_tree_ref>` with it. +While :doc:`RipsComplex <rips_complex_ref>` is convenient, for instance to build a simplicial complex in one line + +.. testcode:: -Finally, it is asked to display information about the Rips complex. + import gudhi + points = [[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]] + cplx = gudhi.RipsComplex(points=points, max_edge_length=12.0).create_simplex_tree(max_dimension=2) +you can achieve the same result without this class for more flexibility .. testcode:: - import gudhi - off_file = gudhi.__root_source_dir__ + '/data/points/alphacomplexdoc.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) - result_str = 'Rips 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(): - print(fmt % tuple(filtered_value)) + import gudhi + from scipy.spatial.distance import cdist + points = [[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]] + distance_matrix = cdist(points, points) + cplx = gudhi.SimplexTree.create_from_array(distance_matrix, max_filtration=12.0) + cplx.expansion(2) -the program output is: +or -.. testoutput:: +.. testcode:: + + import gudhi + from scipy.spatial import cKDTree + points = [[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]] + tree = cKDTree(points) + edges = tree.sparse_distance_matrix(tree, max_distance=12.0, output_type="coo_matrix") + cplx = gudhi.SimplexTree() + cplx.insert_edges_from_coo_matrix(edges) + cplx.expansion(2) - Rips complex is of dimension 1 - 18 simplices - 7 vertices. - [0] -> 0.00 - [1] -> 0.00 - [2] -> 0.00 - [3] -> 0.00 - [4] -> 0.00 - [5] -> 0.00 - [6] -> 0.00 - [2, 3] -> 5.00 - [4, 5] -> 5.39 - [0, 2] -> 5.83 - [0, 1] -> 6.08 - [1, 3] -> 6.32 - [1, 2] -> 6.71 - [5, 6] -> 7.28 - [2, 4] -> 8.94 - [0, 3] -> 9.43 - [4, 6] -> 9.49 - [3, 6] -> 11.00 + +This way, you can easily add a call to :func:`~gudhi.SimplexTree.collapse_edges` before the expansion, +use a different metric to compute the matrix, or other variations. Distance matrix --------------- @@ -223,54 +217,7 @@ until dimension 1 - one skeleton graph in other words), the output is: [4, 6] -> 9.49 [3, 6] -> 11.00 -Example from csv file -^^^^^^^^^^^^^^^^^^^^^ - -This example builds the :doc:`RipsComplex <rips_complex_ref>` from the given -distance matrix in a csv file, and max_edge_length value. -Then it creates a :doc:`SimplexTree <simplex_tree_ref>` with it. - -Finally, it is asked to display information about the Rips complex. - - -.. testcode:: - - import gudhi - distance_matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file=gudhi.__root_source_dir__ + \ - '/data/distance_matrix/full_square_distance_matrix.csv') - rips_complex = gudhi.RipsComplex(distance_matrix=distance_matrix, max_edge_length=12.0) - simplex_tree = rips_complex.create_simplex_tree(max_dimension=1) - result_str = 'Rips 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(): - print(fmt % tuple(filtered_value)) - -the program output is: - -.. testoutput:: - - Rips complex is of dimension 1 - 18 simplices - 7 vertices. - [0] -> 0.00 - [1] -> 0.00 - [2] -> 0.00 - [3] -> 0.00 - [4] -> 0.00 - [5] -> 0.00 - [6] -> 0.00 - [2, 3] -> 5.00 - [4, 5] -> 5.39 - [0, 2] -> 5.83 - [0, 1] -> 6.08 - [1, 3] -> 6.32 - [1, 2] -> 6.71 - [5, 6] -> 7.28 - [2, 4] -> 8.94 - [0, 3] -> 9.43 - [4, 6] -> 9.49 - [3, 6] -> 11.00 +In case this lower triangular matrix is stored in a CSV file, like `data/distance_matrix/full_square_distance_matrix.csv` in the Gudhi distribution, you can read it with :func:`~gudhi.read_lower_triangular_matrix_from_csv_file`. Correlation matrix ------------------ diff --git a/src/python/gudhi/hera/bottleneck.cc b/src/python/gudhi/hera/bottleneck.cc index 0cb562ce..ec461f7c 100644 --- a/src/python/gudhi/hera/bottleneck.cc +++ b/src/python/gudhi/hera/bottleneck.cc @@ -16,7 +16,7 @@ using py::ssize_t; #endif -#include <bottleneck.h> // Hera +#include <hera/bottleneck.h> // Hera double bottleneck_distance(Dgm d1, Dgm d2, double delta) { diff --git a/src/python/gudhi/hera/wasserstein.cc b/src/python/gudhi/hera/wasserstein.cc index fa0cf8aa..3516352e 100644 --- a/src/python/gudhi/hera/wasserstein.cc +++ b/src/python/gudhi/hera/wasserstein.cc @@ -8,10 +8,16 @@ * - YYYY/MM Author: Description of the modification */ -#include <wasserstein.h> // Hera - #include <pybind11_diagram_utils.h> +#ifdef _MSC_VER +// https://github.com/grey-narn/hera/issues/3 +// ssize_t is a non-standard type (well, posix) +using py::ssize_t; +#endif + +#include <hera/wasserstein.h> // Hera + double wasserstein_distance( Dgm d1, Dgm d2, double wasserstein_power, double internal_p, diff --git a/src/python/gudhi/off_reader.pyx b/src/python/gudhi/off_utils.pyx index a3200704..9276c7b0 100644 --- a/src/python/gudhi/off_reader.pyx +++ b/src/python/gudhi/off_utils.pyx @@ -13,8 +13,10 @@ from __future__ import print_function from cython cimport numeric from libcpp.vector cimport vector from libcpp.string cimport string +cimport cython import errno import os +import numpy as np __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" @@ -24,7 +26,7 @@ cdef extern from "Off_reader_interface.h" namespace "Gudhi": vector[vector[double]] read_points_from_OFF_file(string off_file) def read_points_from_off_file(off_file=''): - """Read points from OFF file. + """Read points from an `OFF file <fileformats.html#off-file-format>`_. :param off_file: An OFF file style name. :type off_file: string @@ -39,3 +41,22 @@ def read_points_from_off_file(off_file=''): raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), off_file) +@cython.embedsignature(True) +def write_points_to_off_file(fname, points): + """Write points to an `OFF file <fileformats.html#off-file-format>`_. + + A simple wrapper for `numpy.savetxt`. + + :param fname: Name of the OFF file. + :type fname: str or file handle + :param points: Point coordinates. + :type points: numpy array of shape (n, dim) + """ + points = np.array(points, copy=False) + assert len(points.shape) == 2 + dim = points.shape[1] + if dim == 3: + head = 'OFF\n{} 0 0'.format(points.shape[0]) + else: + head = 'nOFF\n{} {} 0 0'.format(dim, points.shape[0]) + np.savetxt(fname, points, header=head, comments='') diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py index 21275cdd..e438aa66 100644 --- a/src/python/gudhi/persistence_graphical_tools.py +++ b/src/python/gudhi/persistence_graphical_tools.py @@ -194,19 +194,21 @@ def plot_persistence_barcode( y=[(death - birth) if death != float("inf") else (infinity - birth) for (dim,(birth,death)) in persistence] c=[colormap[dim] for (dim,(birth,death)) in persistence] - axes.barh(list(reversed(range(len(x)))), y, height=0.8, left=x, alpha=alpha, color=c, linewidth=0) + axes.barh(range(len(x)), y, left=x, alpha=alpha, color=c, linewidth=0) if legend: - dimensions = list(set(item[0] for item in persistence)) + dimensions = set(item[0] for item in persistence) axes.legend( handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions], loc="lower right", ) axes.set_title("Persistence barcode", fontsize=fontsize) + axes.set_yticks([]) + axes.invert_yaxis() # Ends plot on infinity value and starts a little bit before min_birth if len(x) != 0: - axes.axis([axis_start, infinity, 0, len(x)]) + axes.set_xlim((axis_start, infinity)) return axes except ImportError as import_error: diff --git a/src/python/gudhi/point_cloud/knn.py b/src/python/gudhi/point_cloud/knn.py index de5844f9..7dc83817 100644 --- a/src/python/gudhi/point_cloud/knn.py +++ b/src/python/gudhi/point_cloud/knn.py @@ -314,7 +314,9 @@ class KNearestNeighbors: return None if self.params["implementation"] == "ckdtree": - qargs = {key: val for key, val in self.params.items() if key in {"p", "eps", "n_jobs"}} + qargs = {key: val for key, val in self.params.items() if key in {"p", "eps"}} + # SciPy renamed n_jobs to workers + qargs["workers"] = self.params.get("workers") or self.params.get("n_jobs") or 1 distances, neighbors = self.kdtree.query(X, k=self.k, **qargs) if k == 1: # SciPy decided to squeeze the last dimension for k=1 diff --git a/src/python/gudhi/representations/vector_methods.py b/src/python/gudhi/representations/vector_methods.py index 69ff5e1e..ce74aee5 100644 --- a/src/python/gudhi/representations/vector_methods.py +++ b/src/python/gudhi/representations/vector_methods.py @@ -13,8 +13,13 @@ 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 +try: + # New location since 1.0 + from sklearn.metrics import DistanceMetric +except ImportError: + # Will be removed in 1.3 + from sklearn.neighbors import DistanceMetric from .preprocessing import DiagramScaler, BirthPersistenceTransform @@ -85,7 +90,7 @@ class PersistenceImage(BaseEstimator, TransformerMixin): Xfit.append(image.flatten()[np.newaxis,:]) - Xfit = np.concatenate(Xfit,0) + Xfit = np.concatenate(Xfit, 0) return Xfit @@ -101,7 +106,7 @@ class PersistenceImage(BaseEstimator, TransformerMixin): """ return self.fit_transform([diag])[0,:] -def _automatic_sample_range(sample_range, X, y): +def _automatic_sample_range(sample_range, X): """ Compute and returns sample range from the persistence diagrams if one of the sample_range values is numpy.nan. @@ -114,7 +119,7 @@ def _automatic_sample_range(sample_range, X, y): 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) + pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(X) [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) @@ -123,11 +128,35 @@ def _automatic_sample_range(sample_range, X, y): pass return sample_range + +def _trim_endpoints(x, are_endpoints_nan): + if are_endpoints_nan[0]: + x = x[1:] + if are_endpoints_nan[1]: + x = x[:-1] + return x + + +def _grid_from_sample_range(self, X): + sample_range = np.array(self.sample_range) + self.nan_in_range = np.isnan(sample_range) + self.new_resolution = self.resolution + if not self.keep_endpoints: + self.new_resolution += self.nan_in_range.sum() + self.sample_range_fixed = _automatic_sample_range(sample_range, X) + self.grid_ = np.linspace(self.sample_range_fixed[0], self.sample_range_fixed[1], self.new_resolution) + if not self.keep_endpoints: + self.grid_ = _trim_endpoints(self.grid_, self.nan_in_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. + + Attributes: + grid_ (1d array): The grid on which the landscapes are computed. """ - def __init__(self, num_landscapes=5, resolution=100, sample_range=[np.nan, np.nan]): + def __init__(self, num_landscapes=5, resolution=100, sample_range=[np.nan, np.nan], *, keep_endpoints=False): """ Constructor for the Landscape class. @@ -135,10 +164,10 @@ class Landscape(BaseEstimator, TransformerMixin): num_landscapes (int): number of piecewise-linear functions to output (default 5). resolution (int): number of sample for all piecewise-linear functions (default 100). sample_range ([double, double]): minimum and maximum of all piecewise-linear function domains, 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. + keep_endpoints (bool): when computing `sample_range`, use the exact extremities (where the value is always 0). This is mostly useful for plotting, the default is to use a slightly smaller range. """ self.num_landscapes, self.resolution, self.sample_range = num_landscapes, resolution, sample_range - self.nan_in_range = np.isnan(np.array(self.sample_range)) - self.new_resolution = self.resolution + self.nan_in_range.sum() + self.keep_endpoints = keep_endpoints def fit(self, X, y=None): """ @@ -148,7 +177,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). """ - self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y) + _grid_from_sample_range(self, X) return self def transform(self, X): @@ -161,53 +190,26 @@ class Landscape(BaseEstimator, TransformerMixin): Returns: numpy array with shape (number of diagrams) x (number of samples = **num_landscapes** x **resolution**): output persistence landscapes. """ - num_diag, Xfit = len(X), [] - x_values = np.linspace(self.sample_range[0], self.sample_range[1], self.new_resolution) - step_x = x_values[1] - x_values[0] - - for i in range(num_diag): - - diagram, num_pts_in_diag = X[i], X[i].shape[0] - - ls = np.zeros([self.num_landscapes, self.new_resolution]) - events = [] - for j in range(self.new_resolution): - events.append([]) + Xfit = [] + x_values = self.grid_ + for diag in X: + midpoints, heights = (diag[:, 0] + diag[:, 1]) / 2., (diag[:, 1] - diag[:, 0]) / 2. + tent_functions = np.maximum(heights[None, :] - np.abs(x_values[:, None] - midpoints[None, :]), 0) + n_points = diag.shape[0] + # Complete the array with zeros to get the right number of landscapes + if self.num_landscapes > n_points: + tent_functions = np.concatenate( + [tent_functions, np.zeros((tent_functions.shape[0], self.num_landscapes-n_points))], + axis=1 + ) + tent_functions.partition(tent_functions.shape[1]-self.num_landscapes, axis=1) + landscapes = np.sort(tent_functions[:, -self.num_landscapes:], axis=1)[:, ::-1].T - 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.new_resolution) - mid_idx = np.clip(np.ceil((0.5*(py+px) - self.sample_range[0]) / step_x).astype(int), 0, self.new_resolution) - max_idx = np.clip(np.ceil((py - self.sample_range[0]) / step_x).astype(int), 0, self.new_resolution) - - if min_idx < self.new_resolution and max_idx > 0: - - landscape_value = self.sample_range[0] + min_idx * step_x - px - for k in range(min_idx, mid_idx): - events[k].append(landscape_value) - landscape_value += step_x - - landscape_value = py - self.sample_range[0] - mid_idx * step_x - for k in range(mid_idx, max_idx): - events[k].append(landscape_value) - landscape_value -= step_x + landscapes = np.sqrt(2) * np.ravel(landscapes) + Xfit.append(landscapes) - for j in range(self.new_resolution): - events[j].sort(reverse=True) - for k in range( min(self.num_landscapes, len(events[j])) ): - ls[k,j] = events[j][k] - - if self.nan_in_range[0]: - ls = ls[:,1:] - if self.nan_in_range[1]: - ls = ls[:,:-1] - ls = np.sqrt(2)*np.reshape(ls,[1,-1]) - Xfit.append(ls) - - Xfit = np.concatenate(Xfit,0) - - return Xfit + return np.stack(Xfit, axis=0) def __call__(self, diag): """ @@ -219,13 +221,16 @@ class Landscape(BaseEstimator, TransformerMixin): Returns: numpy array with shape (number of samples = **num_landscapes** x **resolution**): output persistence landscape. """ - return self.fit_transform([diag])[0,:] + return self.fit_transform([diag])[0, :] class Silhouette(BaseEstimator, TransformerMixin): """ This is a class for computing persistence silhouettes from a list of persistence diagrams. A persistence silhouette is computed by taking a weighted average of the collection of 1D piecewise-linear functions given by the persistence landscapes, and then by evenly sampling this average on a given range. Finally, the corresponding vector of samples is returned. See https://arxiv.org/abs/1312.0308 for more details. + + Attributes: + grid_ (1d array): The grid on which the silhouette is computed. """ - def __init__(self, weight=lambda x: 1, resolution=100, sample_range=[np.nan, np.nan]): + def __init__(self, weight=lambda x: 1, resolution=100, sample_range=[np.nan, np.nan], *, keep_endpoints=False): """ Constructor for the Silhouette class. @@ -233,8 +238,10 @@ class Silhouette(BaseEstimator, TransformerMixin): weight (function): weight function for the persistence diagram points (default constant function, ie lambda x: 1). This function must be defined on 2D points, ie on lists or numpy arrays of the form [p_x,p_y]. resolution (int): number of samples for the weighted average (default 100). sample_range ([double, double]): minimum and maximum for the weighted average 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. + keep_endpoints (bool): when computing `sample_range`, use the exact extremities (where the value is always 0). This is mostly useful for plotting, the default is to use a slightly smaller range. """ self.weight, self.resolution, self.sample_range = weight, resolution, sample_range + self.keep_endpoints = keep_endpoints def fit(self, X, y=None): """ @@ -244,7 +251,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). """ - self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y) + _grid_from_sample_range(self, X) return self def transform(self, X): @@ -257,44 +264,19 @@ class Silhouette(BaseEstimator, TransformerMixin): Returns: numpy array with shape (number of diagrams) x (**resolution**): output persistence silhouettes. """ - 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): - - diagram, num_pts_in_diag = X[i], X[i].shape[0] + Xfit = [] + x_values = self.grid_ - sh, weights = np.zeros(self.resolution), np.zeros(num_pts_in_diag) - for j in range(num_pts_in_diag): - weights[j] = self.weight(diagram[j,:]) + for diag in X: + midpoints, heights = (diag[:, 0] + diag[:, 1]) / 2., (diag[:, 1] - diag[:, 0]) / 2. + weights = np.array([self.weight(pt) for pt in diag]) total_weight = np.sum(weights) - for j in range(num_pts_in_diag): - - [px,py] = diagram[j,:2] - weight = weights[j] / total_weight - min_idx = np.clip(np.ceil((px - self.sample_range[0]) / step_x).astype(int), 0, self.resolution) - mid_idx = np.clip(np.ceil((0.5*(py+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) - - if min_idx < self.resolution and max_idx > 0: + tent_functions = np.maximum(heights[None, :] - np.abs(x_values[:, None] - midpoints[None, :]), 0) + silhouette = np.sum(weights[None, :] / total_weight * tent_functions, axis=1) + Xfit.append(silhouette * np.sqrt(2)) - silhouette_value = self.sample_range[0] + min_idx * step_x - px - for k in range(min_idx, mid_idx): - sh[k] += weight * silhouette_value - silhouette_value += step_x - - silhouette_value = py - self.sample_range[0] - mid_idx * step_x - for k in range(mid_idx, max_idx): - sh[k] += weight * silhouette_value - silhouette_value -= step_x - - Xfit.append(np.reshape(np.sqrt(2) * sh, [1,-1])) - - Xfit = np.concatenate(Xfit, 0) - - return Xfit + return np.stack(Xfit, axis=0) def __call__(self, diag): """ @@ -312,36 +294,39 @@ class Silhouette(BaseEstimator, TransformerMixin): class BettiCurve(BaseEstimator, TransformerMixin): """ 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], predefined_grid=None): - """ - Constructor for the BettiCurve class. + 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: - 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. + >>> bc = BettiCurve(predefined_grid=xs) # doctest: +SKIP + >>> result = bc(pd) # doctest: +SKIP - 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. + and - 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: + >>> 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 - >>> bc = BettiCurve(predefined_grid=xs) # doctest: +SKIP - >>> result = bc(pd) # doctest: +SKIP + are the same. - and + 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 and resolution is None, the grid is fitted to capture all filtration values at which the Betti numbers change. + """ - >>> 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 + def __init__(self, resolution=100, sample_range=[np.nan, np.nan], predefined_grid=None, *, keep_endpoints=False): + """ + Constructor for the BettiCurve class. - are the same. + Parameters: + resolution (int): number of samples for the piecewise-constant function (default 100), or None for the exact curve. + 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. + keep_endpoints (bool): when computing `sample_range` (fixed `resolution`, no `predefined_grid`), use the exact extremities. This is mostly useful for plotting, the default is to use a slightly smaller range. """ if (predefined_grid is not None) and (not isinstance(predefined_grid, np.ndarray)): @@ -350,6 +335,7 @@ class BettiCurve(BaseEstimator, TransformerMixin): self.predefined_grid = predefined_grid self.resolution = resolution self.sample_range = sample_range + self.keep_endpoints = keep_endpoints def is_fitted(self): return hasattr(self, "grid_") @@ -368,8 +354,7 @@ class BettiCurve(BaseEstimator, TransformerMixin): 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) + _grid_from_sample_range(self, X) else: self.grid_ = self.predefined_grid # Get the predefined grid from user @@ -468,8 +453,11 @@ class BettiCurve(BaseEstimator, TransformerMixin): class Entropy(BaseEstimator, TransformerMixin): """ This is a class for computing persistence entropy. Persistence entropy is a statistic for persistence diagrams inspired from Shannon entropy. This statistic can also be used to compute a feature vector, called the entropy summary function. See https://arxiv.org/pdf/1803.08304.pdf for more details. Note that a previous implementation was contributed by Manuel Soriano-Trigueros. + + Attributes: + grid_ (1d array): In vector mode, the grid on which the entropy summary function is computed. """ - def __init__(self, mode="scalar", normalized=True, resolution=100, sample_range=[np.nan, np.nan]): + def __init__(self, mode="scalar", normalized=True, resolution=100, sample_range=[np.nan, np.nan], *, keep_endpoints=False): """ Constructor for the Entropy class. @@ -478,8 +466,10 @@ class Entropy(BaseEstimator, TransformerMixin): normalized (bool): whether to normalize the entropy summary function (default True). Used only if **mode** = "vector". resolution (int): number of sample for the entropy summary function (default 100). Used only if **mode** = "vector". sample_range ([double, double]): minimum and maximum of the entropy summary 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. Used only if **mode** = "vector". + keep_endpoints (bool): when computing `sample_range`, use the exact extremities. This is mostly useful for plotting, the default is to use a slightly smaller range. """ self.mode, self.normalized, self.resolution, self.sample_range = mode, normalized, resolution, sample_range + self.keep_endpoints = keep_endpoints def fit(self, X, y=None): """ @@ -489,7 +479,9 @@ class Entropy(BaseEstimator, TransformerMixin): X (list of n x 2 numpy arrays): input persistence diagrams. y (n x 1 array): persistence diagram labels (unused). """ - self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y) + if self.mode == "vector": + _grid_from_sample_range(self, X) + self.step_ = self.grid_[1] - self.grid_[0] return self def transform(self, X): @@ -503,8 +495,6 @@ class Entropy(BaseEstimator, TransformerMixin): numpy array with shape (number of diagrams) x (1 if **mode** = "scalar" else **resolution**): output entropy. """ 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] new_X = BirthPersistenceTransform().fit_transform(X) for i in range(num_diag): @@ -519,8 +509,8 @@ class Entropy(BaseEstimator, TransformerMixin): ent = np.zeros(self.resolution) for j in range(num_pts_in_diag): [px,py] = orig_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) + min_idx = np.clip(np.ceil((px - self.sample_range_fixed[0]) / self.step_).astype(int), 0, self.resolution) + max_idx = np.clip(np.ceil((py - self.sample_range_fixed[0]) / self.step_).astype(int), 0, self.resolution) ent[min_idx:max_idx]-=p[j]*np.log(p[j]) if self.normalized: ent = ent / np.linalg.norm(ent, ord=1) @@ -720,17 +710,17 @@ class Atol(BaseEstimator, TransformerMixin): >>> 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.fit(X=[a, b, c]).centers + array([[ 2.6 , 2.8 , -0.4 ], + [ 2. , 0.66666667, 3.33333333]]) >>> atol_vectoriser(a) - >>> # array([1.18168665, 0.42375966]) # doctest: +SKIP + array([0.42375966, 1.18168665]) >>> 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]]) + array([1.25157463, 0.02062512]) + >>> atol_vectoriser.transform(X=[a, b, c]) + array([[0.42375966, 1.18168665], + [1.06330156, 0.29861028], + [1.25157463, 0.02062512]]) """ # 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"): @@ -781,6 +771,8 @@ class Atol(BaseEstimator, TransformerMixin): measures_concat = np.concatenate(X) self.quantiser.fit(X=measures_concat, sample_weight=sample_weight) self.centers = self.quantiser.cluster_centers_ + # Hack, but some people are unhappy if the order depends on the version of sklearn + self.centers = self.centers[np.lexsort(self.centers.T)] if self.quantiser.n_clusters == 1: dist_centers = pairwise.pairwise_distances(measures_concat) np.fill_diagonal(dist_centers, 0) diff --git a/src/python/gudhi/rips_complex.pyx b/src/python/gudhi/rips_complex.pyx index c3470292..d748f91e 100644 --- a/src/python/gudhi/rips_complex.pyx +++ b/src/python/gudhi/rips_complex.pyx @@ -41,31 +41,30 @@ cdef class RipsComplex: cdef Rips_complex_interface thisref # Fake constructor that does nothing but documenting the constructor - def __init__(self, points=None, distance_matrix=None, + def __init__(self, *, points=None, distance_matrix=None, max_edge_length=float('inf'), sparse=None): """RipsComplex constructor. - :param max_edge_length: Rips value. - :type max_edge_length: float - :param points: A list of points in d-Dimension. - :type points: list of list of float + :type points: List[List[float]] Or :param distance_matrix: A distance matrix (full square or lower triangular). - :type points: list of list of float + :type distance_matrix: List[List[float]] And in both cases + :param max_edge_length: Rips value. + :type max_edge_length: float :param sparse: If this is not None, it switches to building a sparse Rips and represents the approximation parameter epsilon. :type sparse: float """ # The real cython constructor - def __cinit__(self, points=None, distance_matrix=None, + def __cinit__(self, *, points=None, distance_matrix=None, max_edge_length=float('inf'), sparse=None): if sparse is not None: if distance_matrix is not None: diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 5642f82d..5309c6fa 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -56,6 +56,8 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": int upper_bound_dimension() nogil bool find_simplex(vector[int] simplex) nogil bool insert(vector[int] simplex, double filtration) nogil + void insert_matrix(double* filtrations, int n, int stride0, int stride1, double max_filtration) nogil except + + void insert_batch_vertices(vector[int] v, double f) nogil except + vector[pair[vector[int], double]] get_star(vector[int] simplex) nogil vector[pair[vector[int], double]] get_cofaces(vector[int] simplex, int dimension) nogil void expansion(int max_dim) nogil except + diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 05bfe22e..4cf176f5 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -8,14 +8,24 @@ # - YYYY/MM Author: Description of the modification from cython.operator import dereference, preincrement -from libc.stdint cimport intptr_t +from libc.stdint cimport intptr_t, int32_t, int64_t import numpy as np cimport gudhi.simplex_tree +cimport cython +from numpy.math cimport INFINITY __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" __license__ = "MIT" +ctypedef fused some_int: + int32_t + int64_t + +ctypedef fused some_float: + float + double + cdef bool callback(vector[int] simplex, void *blocker_func): return (<object>blocker_func)(simplex) @@ -228,6 +238,91 @@ cdef class SimplexTree: """ return self.get_ptr().insert(simplex, <double>filtration) + @staticmethod + @cython.boundscheck(False) + def create_from_array(filtrations, double max_filtration=INFINITY): + """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:: + + filtrations[np.diag_indices_from(filtrations)] = filtrations.min(axis=1) + + or:: + + diag = filtrations.diagonal() + filtrations = np.fmax(np.fmax(filtrations, diag[:, None]), diag[None, :]) + + :param filtrations: the filtration values of the vertices and edges to insert. The matrix is assumed to be symmetric. + :type filtrations: numpy.ndarray of shape (n,n) + :param max_filtration: only insert vertices and edges with filtration values no larger than max_filtration + :type max_filtration: float + :returns: the new complex + :rtype: SimplexTree + """ + # TODO: document which half of the matrix is actually read? + filtrations = np.asanyarray(filtrations, dtype=float) + cdef double[:,:] F = filtrations + ret = SimplexTree() + cdef int n = F.shape[0] + assert n == F.shape[1], 'create_from_array() expects a square array' + with nogil: + ret.get_ptr().insert_matrix(&F[0,0], n, F.strides[0], F.strides[1], max_filtration) + return ret + + def insert_edges_from_coo_matrix(self, edges): + """Inserts edges given by a sparse matrix in `COOrdinate format + <https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html>`_. + If an edge is repeated, the smallest filtration value is used. Missing entries are not inserted. + Diagonal entries are currently interpreted as vertices, although we do not guarantee this behavior + in the future, and this is only useful if you want to insert vertices with a smaller filtration value + than the smallest edge containing it, since vertices are implicitly inserted together with the edges. + + :param edges: the edges to insert and their filtration values. + :type edges: scipy.sparse.coo_matrix of shape (n,n) + + .. seealso:: :func:`insert_batch` + """ + # Without this, it could be slow if we end up inserting vertices in a bad order (flat_map). + self.get_ptr().insert_batch_vertices(np.unique(np.stack((edges.row, edges.col))), INFINITY) + # TODO: optimize this? + for edge in zip(edges.row, edges.col, edges.data): + self.get_ptr().insert((edge[0], edge[1]), edge[2]) + + @cython.boundscheck(False) + @cython.wraparound(False) + def insert_batch(self, some_int[:,:] vertex_array, some_float[:] filtrations): + """Inserts k-simplices given by a sparse array in a format similar + to `torch.sparse <https://pytorch.org/docs/stable/sparse.html>`_. + The n-th simplex has vertices `vertex_array[0,n]`, ..., + `vertex_array[k,n]` and filtration value `filtrations[n]`. + If a simplex is repeated, the smallest filtration value is used. + Simplices with a repeated vertex are currently interpreted as lower + dimensional simplices, but we do not guarantee this behavior in the + future. Any time a simplex is inserted, its faces are inserted as well + if needed to preserve a simplicial complex. + + :param vertex_array: the k-simplices to insert. + :type vertex_array: numpy.array of shape (k+1,n) + :param filtrations: the filtration values. + :type filtrations: numpy.array of shape (n,) + """ + cdef vector[int] vertices = np.unique(vertex_array) + cdef Py_ssize_t k = vertex_array.shape[0] + cdef Py_ssize_t n = vertex_array.shape[1] + assert filtrations.shape[0] == n, 'inconsistent sizes for vertex_array and filtrations' + cdef Py_ssize_t i + cdef Py_ssize_t j + cdef vector[int] v + with nogil: + # Without this, it could be slow if we end up inserting vertices in a bad order (flat_map). + # NaN currently does the wrong thing + self.get_ptr().insert_batch_vertices(vertices, INFINITY) + for i in range(n): + for j in range(k): + v.push_back(vertex_array[j, i]) + self.get_ptr().insert(v, filtrations[i]) + v.clear() + def get_simplices(self): """This function returns a generator with simplices and their given filtration values. @@ -376,7 +471,7 @@ cdef class SimplexTree: """ return self.get_ptr().prune_above_filtration(filtration) - def expansion(self, max_dim): + def expansion(self, max_dimension): """Expands the simplex tree containing only its one skeleton until dimension max_dim. @@ -390,10 +485,10 @@ cdef class SimplexTree: 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 + :param max_dimension: The maximal dimension. + :type max_dimension: int """ - cdef int maxdim = max_dim + cdef int maxdim = max_dimension with nogil: self.get_ptr().expansion(maxdim) diff --git a/src/python/gudhi/tensorflow/cubical_layer.py b/src/python/gudhi/tensorflow/cubical_layer.py index 3304e719..5df2c370 100644 --- a/src/python/gudhi/tensorflow/cubical_layer.py +++ b/src/python/gudhi/tensorflow/cubical_layer.py @@ -18,7 +18,7 @@ def _Cubical(Xflat, Xdim, dimensions, homology_coeff_field): cc = CubicalComplex(dimensions=Xdim[::-1], top_dimensional_cells=Xflat) cc.compute_persistence(homology_coeff_field=homology_coeff_field) - # Retrieve and ouput image indices/pixels corresponding to positive and negative simplices + # Retrieve and output image indices/pixels corresponding to positive and negative simplices cof_pp = cc.cofaces_of_persistence_pairs() L_cofs = [] diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index 3d20aa8f..41eb72c1 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -106,7 +106,7 @@ class Exact_alpha_complex_dD final : 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 { + virtual std::size_t num_vertices() const override { return alpha_complex_.num_vertices(); } @@ -141,7 +141,7 @@ class Inexact_alpha_complex_dD final : public Abstract_alpha_complex { return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, false, default_filtration_value); } - virtual std::size_t num_vertices() const { + virtual std::size_t num_vertices() const override { return alpha_complex_.num_vertices(); } diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 3848c5ad..0317ea39 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -40,6 +40,8 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { 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; typedef bool (*blocker_func_t)(Simplex simplex, void *user_data); public: @@ -62,6 +64,30 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { return (result.second); } + void insert_matrix(double* filtrations, int n, int stride0, int stride1, double max_filtration) { + // We could delegate to insert_graph, but wrapping the matrix in a graph interface is too much work, + // and this is a bit more efficient. + auto& rm = this->root()->members_; + for(int i=0; i<n; ++i) { + char* p = reinterpret_cast<char*>(filtrations) + i * stride0; + double fv = *reinterpret_cast<double*>(p + i * stride1); + if(fv > max_filtration) continue; + auto sh = rm.emplace_hint(rm.end(), i, Node(this->root(), fv)); + Siblings* children = nullptr; + // Should we make a first pass to count the number of edges so we can reserve the right space? + for(int j=i+1; j<n; ++j) { + double fe = *reinterpret_cast<double*>(p + j * stride1); + if(fe > max_filtration) continue; + if(!children) { + children = new Siblings(this->root(), i); + sh->second.assign_children(children); + } + children->members().emplace_hint(children->members().end(), j, Node(children, fe)); + } + } + + } + // Do not interface this function, only used in alpha complex interface for complex creation bool insert_simplex(const Simplex& simplex, Filtration_value filtration = 0) { Insertion_result result = Base::insert_simplex(simplex, filtration); diff --git a/src/python/setup.py.in b/src/python/setup.py.in index 2c67c2c5..6eb0db42 100644 --- a/src/python/setup.py.in +++ b/src/python/setup.py.in @@ -48,10 +48,6 @@ 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)] - if module == 'hera/wasserstein': - my_include_dirs = ['@HERA_WASSERSTEIN_INCLUDE_DIR@'] + my_include_dirs - elif module == 'hera/bottleneck': - my_include_dirs = ['@HERA_BOTTLENECK_INCLUDE_DIR@'] + my_include_dirs ext_modules.append(Extension( 'gudhi.' + module.replace('/', '.'), sources = [source_dir + module + '.cc'], diff --git a/src/python/test/test_off.py b/src/python/test/test_off.py new file mode 100644 index 00000000..aea1941b --- /dev/null +++ b/src/python/test/test_off.py @@ -0,0 +1,21 @@ +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + Author(s): Marc Glisse + + Copyright (C) 2022 Inria + + Modification(s): + - YYYY/MM Author: Description of the modification +""" + +import gudhi as gd +import numpy as np +import pytest + + +def test_off_rw(): + for dim in range(2, 6): + X = np.random.rand(123, dim) + gd.write_points_to_off_file("rand.off", X) + Y = gd.read_points_from_off_file("rand.off") + assert Y == pytest.approx(X) diff --git a/src/python/test/test_persistence_graphical_tools.py b/src/python/test/test_persistence_graphical_tools.py index c19836b7..0e2ac3f8 100644 --- a/src/python/test/test_persistence_graphical_tools.py +++ b/src/python/test/test_persistence_graphical_tools.py @@ -12,6 +12,7 @@ import gudhi as gd import numpy as np import matplotlib as plt import pytest +import warnings def test_array_handler(): @@ -71,13 +72,13 @@ def test_limit_to_max_intervals(): (0, (0.0, 0.106382)), ] # check no warnings if max_intervals equals to the diagrams number - with pytest.warns(None) as record: + with warnings.catch_warnings(): + warnings.simplefilter("error") truncated_diags = gd.persistence_graphical_tools._limit_to_max_intervals( diags, 10, key=lambda life_time: life_time[1][1] - life_time[1][0] ) # check diagrams are not sorted assert truncated_diags == diags - assert len(record) == 0 # check warning if max_intervals lower than the diagrams number with pytest.warns(UserWarning) as record: diff --git a/src/python/test/test_representations.py b/src/python/test/test_representations.py index 4a455bb6..f4ffbdc1 100755 --- a/src/python/test/test_representations.py +++ b/src/python/test/test_representations.py @@ -161,7 +161,7 @@ def test_entropy_miscalculation(): return -np.dot(l, np.log(l)) sce = Entropy(mode="scalar") assert [[pe(diag_ex)]] == sce.fit_transform([diag_ex]) - sce = Entropy(mode="vector", resolution=4, normalized=False) + sce = Entropy(mode="vector", resolution=4, normalized=False, keep_endpoints=True) pef = [-1/4*np.log(1/4)-1/4*np.log(1/4)-1/2*np.log(1/2), -1/4*np.log(1/4)-1/4*np.log(1/4)-1/2*np.log(1/2), -1/2*np.log(1/2), @@ -170,7 +170,7 @@ def test_entropy_miscalculation(): sce = Entropy(mode="vector", resolution=4, normalized=True) pefN = (sce.fit_transform([diag_ex]))[0] area = np.linalg.norm(pefN, ord=1) - assert area==1 + assert area==pytest.approx(1) def test_kernel_empty_diagrams(): empty_diag = np.empty(shape = [0, 2]) @@ -187,3 +187,83 @@ def test_kernel_empty_diagrams(): # 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) + +def test_silhouette_permutation_invariance(): + dgm = _n_diags(1)[0] + dgm_permuted = dgm[np.random.permutation(dgm.shape[0]).astype(int)] + random_resolution = random.randint(50, 100) * 10 + slt = Silhouette(resolution=random_resolution, weight=pow(2)) + + assert np.all(np.isclose(slt(dgm), slt(dgm_permuted))) + + +def test_silhouette_multiplication_invariance(): + dgm = _n_diags(1)[0] + n_repetitions = np.random.randint(2, high=10) + dgm_augmented = np.repeat(dgm, repeats=n_repetitions, axis=0) + + random_resolution = random.randint(50, 100) * 10 + slt = Silhouette(resolution=random_resolution, weight=pow(2)) + assert np.all(np.isclose(slt(dgm), slt(dgm_augmented))) + + +def test_silhouette_numeric(): + dgm = np.array([[2., 3.], [5., 6.]]) + slt = Silhouette(resolution=9, weight=pow(1), sample_range=[2., 6.]) + #slt.fit([dgm]) + # x_values = array([2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6.]) + + expected_silhouette = np.array([0., 0.5, 0., 0., 0., 0., 0., 0.5, 0.])/np.sqrt(2) + output_silhouette = slt(dgm) + assert np.all(np.isclose(output_silhouette, expected_silhouette)) + + +def test_landscape_small_persistence_invariance(): + dgm = np.array([[2., 6.], [2., 5.], [3., 7.]]) + small_persistence_pts = np.random.rand(10, 2) + small_persistence_pts[:, 1] += small_persistence_pts[:, 0] + small_persistence_pts += np.min(dgm) + dgm_augmented = np.concatenate([dgm, small_persistence_pts], axis=0) + + lds = Landscape(num_landscapes=2, resolution=5) + lds_dgm, lds_dgm_augmented = lds(dgm), lds(dgm_augmented) + + assert np.all(np.isclose(lds_dgm, lds_dgm_augmented)) + + +def test_landscape_numeric(): + dgm = np.array([[2., 6.], [3., 5.]]) + lds_ref = np.array([ + 0., 0.5, 1., 1.5, 2., 1.5, 1., 0.5, 0., # tent of [2, 6] + 0., 0., 0., 0.5, 1., 0.5, 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., + ]) + lds_ref *= np.sqrt(2) + lds = Landscape(num_landscapes=4, resolution=9, sample_range=[2., 6.]) + lds_dgm = lds(dgm) + assert np.all(np.isclose(lds_dgm, lds_ref)) + + +def test_landscape_nan_range(): + dgm = np.array([[2., 6.], [3., 5.]]) + lds = Landscape(num_landscapes=2, resolution=9, sample_range=[np.nan, 6.]) + lds_dgm = lds(dgm) + assert (lds.sample_range_fixed[0] == 2) & (lds.sample_range_fixed[1] == 6) + assert lds.new_resolution == 10 + +def test_endpoints(): + diags = [ np.array([[2., 3.]]) ] + for vec in [ Landscape(), Silhouette(), BettiCurve(), Entropy(mode="vector") ]: + vec.fit(diags) + assert vec.grid_[0] > 2 and vec.grid_[-1] < 3 + for vec in [ Landscape(keep_endpoints=True), Silhouette(keep_endpoints=True), BettiCurve(keep_endpoints=True), Entropy(mode="vector", keep_endpoints=True)]: + vec.fit(diags) + assert vec.grid_[0] == 2 and vec.grid_[-1] == 3 + vec = BettiCurve(resolution=None) + vec.fit(diags) + assert np.equal(vec.grid_, [-np.inf, 2., 3.]).all() + +def test_get_params(): + for vec in [ Landscape(), Silhouette(), BettiCurve(), Entropy(mode="vector") ]: + vec.get_params() diff --git a/src/python/test/test_simplex_generators.py b/src/python/test/test_simplex_generators.py index 8a9b4844..c567d4c1 100755 --- a/src/python/test/test_simplex_generators.py +++ b/src/python/test/test_simplex_generators.py @@ -14,7 +14,7 @@ import numpy as np def test_flag_generators(): pts = np.array([[0, 0], [0, 1.01], [1, 0], [1.02, 1.03], [100, 0], [100, 3.01], [103, 0], [103.02, 3.03]]) - r = gudhi.RipsComplex(pts, max_edge_length=4) + r = gudhi.RipsComplex(points=pts, max_edge_length=4) st = r.create_simplex_tree(max_dimension=50) st.persistence() g = st.flag_persistence_generators() diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 54bafed5..2ccbfbf5 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -249,6 +249,7 @@ def test_make_filtration_non_decreasing(): assert st.filtration([3, 4]) == 2.0 assert st.filtration([4, 5]) == 2.0 + def test_extend_filtration(): # Inserted simplex: @@ -257,86 +258,87 @@ def test_extend_filtration(): # / \ / # o o # /2\ /3 - # o o - # 1 0 - - st = SimplexTree() - st.insert([0,2]) - st.insert([1,2]) - st.insert([0,3]) - st.insert([2,5]) - st.insert([3,4]) - st.insert([3,5]) - st.assign_filtration([0], 1.) - st.assign_filtration([1], 2.) - st.assign_filtration([2], 3.) - st.assign_filtration([3], 4.) - st.assign_filtration([4], 5.) - st.assign_filtration([5], 6.) - - assert list(st.get_filtration()) == [ - ([0, 2], 0.0), - ([1, 2], 0.0), - ([0, 3], 0.0), - ([3, 4], 0.0), - ([2, 5], 0.0), - ([3, 5], 0.0), - ([0], 1.0), - ([1], 2.0), - ([2], 3.0), - ([3], 4.0), - ([4], 5.0), - ([5], 6.0) + # o o + # 1 0 + + st = SimplexTree() + st.insert([0, 2]) + st.insert([1, 2]) + st.insert([0, 3]) + st.insert([2, 5]) + st.insert([3, 4]) + st.insert([3, 5]) + st.assign_filtration([0], 1.0) + st.assign_filtration([1], 2.0) + st.assign_filtration([2], 3.0) + st.assign_filtration([3], 4.0) + st.assign_filtration([4], 5.0) + st.assign_filtration([5], 6.0) + + assert list(st.get_filtration()) == [ + ([0, 2], 0.0), + ([1, 2], 0.0), + ([0, 3], 0.0), + ([3, 4], 0.0), + ([2, 5], 0.0), + ([3, 5], 0.0), + ([0], 1.0), + ([1], 2.0), + ([2], 3.0), + ([3], 4.0), + ([4], 5.0), + ([5], 6.0), ] - + st.extend_filtration() - - assert list(st.get_filtration()) == [ - ([6], -3.0), - ([0], -2.0), - ([1], -1.8), - ([2], -1.6), - ([0, 2], -1.6), - ([1, 2], -1.6), - ([3], -1.4), - ([0, 3], -1.4), - ([4], -1.2), - ([3, 4], -1.2), - ([5], -1.0), - ([2, 5], -1.0), - ([3, 5], -1.0), - ([5, 6], 1.0), - ([4, 6], 1.2), - ([3, 6], 1.4), + + assert list(st.get_filtration()) == [ + ([6], -3.0), + ([0], -2.0), + ([1], -1.8), + ([2], -1.6), + ([0, 2], -1.6), + ([1, 2], -1.6), + ([3], -1.4), + ([0, 3], -1.4), + ([4], -1.2), + ([3, 4], -1.2), + ([5], -1.0), + ([2, 5], -1.0), + ([3, 5], -1.0), + ([5, 6], 1.0), + ([4, 6], 1.2), + ([3, 6], 1.4), ([3, 4, 6], 1.4), - ([3, 5, 6], 1.4), - ([2, 6], 1.6), - ([2, 5, 6], 1.6), - ([1, 6], 1.8), - ([1, 2, 6], 1.8), - ([0, 6], 2.0), - ([0, 2, 6], 2.0), - ([0, 3, 6], 2.0) + ([3, 5, 6], 1.4), + ([2, 6], 1.6), + ([2, 5, 6], 1.6), + ([1, 6], 1.8), + ([1, 2, 6], 1.8), + ([0, 6], 2.0), + ([0, 2, 6], 2.0), + ([0, 3, 6], 2.0), ] - dgms = st.extended_persistence(min_persistence=-1.) + dgms = st.extended_persistence(min_persistence=-1.0) assert len(dgms) == 4 # Sort by (death-birth) descending - we are only interested in those with the longest life span for idx in range(4): - dgms[idx] = sorted(dgms[idx], key=lambda x:(-abs(x[1][0]-x[1][1]))) + dgms[idx] = sorted(dgms[idx], key=lambda x: (-abs(x[1][0] - x[1][1]))) + + assert dgms[0][0][1][0] == pytest.approx(2.0) + assert dgms[0][0][1][1] == pytest.approx(3.0) + assert dgms[1][0][1][0] == pytest.approx(5.0) + assert dgms[1][0][1][1] == pytest.approx(4.0) + assert dgms[2][0][1][0] == pytest.approx(1.0) + assert dgms[2][0][1][1] == pytest.approx(6.0) + assert dgms[3][0][1][0] == pytest.approx(6.0) + assert dgms[3][0][1][1] == pytest.approx(1.0) - assert dgms[0][0][1][0] == pytest.approx(2.) - assert dgms[0][0][1][1] == pytest.approx(3.) - assert dgms[1][0][1][0] == pytest.approx(5.) - assert dgms[1][0][1][1] == pytest.approx(4.) - assert dgms[2][0][1][0] == pytest.approx(1.) - assert dgms[2][0][1][1] == pytest.approx(6.) - assert dgms[3][0][1][0] == pytest.approx(6.) - assert dgms[3][0][1][1] == pytest.approx(1.) def test_simplices_iterator(): st = SimplexTree() - + assert st.insert([0, 1, 2], filtration=4.0) == True assert st.insert([2, 3, 4], filtration=2.0) == True @@ -346,9 +348,10 @@ 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 @@ -360,31 +363,33 @@ def test_collapse_edges(): st.collapse_edges() assert st.num_simplices() == 9 - assert st.find([0, 2]) == False # [1, 3] would be fine as well + assert st.find([0, 2]) == False # [1, 3] would be fine as well for simplex in st.get_skeleton(0): - assert simplex[1] == 1. + assert simplex[1] == 1.0 + 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 + + assert st.insert([0, 1, 2], 3.0) == True + assert st.insert([0, 3], 2.0) == True + assert st.insert([3, 4, 5], 3.0) == True + assert st.insert([0, 1, 6, 7], 4.0) == True # Guaranteed by construction for simplex in st.get_simplices(): - assert st.filtration(simplex[0]) >= 2. - + assert st.filtration(simplex[0]) >= 2.0 + # 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) + st.reset_filtration(0.0, dimension) for simplex in st.get_skeleton(3): print(simplex) if len(simplex[0]) < (dimension) + 1: - assert st.filtration(simplex[0]) >= 2. + assert st.filtration(simplex[0]) >= 2.0 else: - assert st.filtration(simplex[0]) == 0. + assert st.filtration(simplex[0]) == 0.0 + def test_boundaries_iterator(): st = SimplexTree() @@ -400,16 +405,17 @@ def test_boundaries_iterator(): list(st.get_boundaries([])) with pytest.raises(RuntimeError): - list(st.get_boundaries([0, 4])) # (0, 4) does not exist + list(st.get_boundaries([0, 4])) # (0, 4) does not exist with pytest.raises(RuntimeError): - list(st.get_boundaries([6])) # (6) does not exist + 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 # | \ | \ | \ | \ | # | \ | \ | \ | \ | @@ -418,50 +424,52 @@ def test_persistence_intervals_in_dimension(): # | \ | \ | \ | \ | # 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.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")]])) + assert np.array_equal(H0, np.array([[0.0, float("inf")]])) H1 = st.persistence_intervals_in_dimension(1) - assert np.array_equal(H1, np.array([[ 0., float("inf")], [ 0., float("inf")]])) + assert np.array_equal(H1, np.array([[0.0, float("inf")], [0.0, float("inf")]])) H2 = st.persistence_intervals_in_dimension(2) - assert np.array_equal(H2, np.array([[ 0., float("inf")]])) + assert np.array_equal(H2, np.array([[0.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.) + st1.insert([1, 2, 3], 4.0) assert st1 != st2 - st2.insert([1,2,3], 4.) + st2.insert([1, 2, 3], 4.0) assert st1 == st2 + def test_simplex_tree_deep_copy(): st = SimplexTree() - st.insert([1, 2, 3], 0.) + st.insert([1, 2, 3], 0.0) # compute persistence only on the original st.compute_persistence() @@ -480,14 +488,15 @@ def test_simplex_tree_deep_copy(): 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.) + st.insert([1, 2, 3], 0.0) # compute persistence only on the original st.compute_persistence() @@ -506,56 +515,132 @@ def test_simplex_tree_deep_copy_constructor(): 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") + 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.0) + 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 + from scipy.spatial import cKDTree + except ImportError: + print("Skipping, no SciPy") + return + + st = SimplexTree() + st.insert([1, 2, 7], 7) + row = np.array([2, 5, 3]) + col = np.array([1, 4, 6]) + dat = np.array([1, 2, 3]) + edges = coo_matrix((dat, (row, col))) + st.insert_edges_from_coo_matrix(edges) + assert list(st.get_filtration()) == [ + ([1], 1.0), + ([2], 1.0), + ([1, 2], 1.0), + ([4], 2.0), + ([5], 2.0), + ([4, 5], 2.0), + ([3], 3.0), + ([6], 3.0), + ([3, 6], 3.0), + ([7], 7.0), + ([1, 7], 7.0), + ([2, 7], 7.0), + ([1, 2, 7], 7.0), + ] + + pts = np.random.rand(100, 2) + tree = cKDTree(pts) + edges = tree.sparse_distance_matrix(tree, max_distance=0.15, output_type="coo_matrix") + st = SimplexTree() + st.insert_edges_from_coo_matrix(edges) + assert 100 < st.num_simplices() < 1000 + + +def test_insert_batch(): + st = SimplexTree() + # vertices + st.insert_batch(np.array([[6, 1, 5]]), np.array([-5.0, 2.0, -3.0])) + # triangles + st.insert_batch(np.array([[2, 10], [5, 0], [6, 11]]), np.array([4.0, 0.0])) + # edges + st.insert_batch(np.array([[1, 5], [2, 5]]), np.array([1.0, 3.0])) + + assert list(st.get_filtration()) == [ + ([6], -5.0), + ([5], -3.0), + ([0], 0.0), + ([10], 0.0), + ([0, 10], 0.0), + ([11], 0.0), + ([0, 11], 0.0), + ([10, 11], 0.0), + ([0, 10, 11], 0.0), + ([1], 1.0), + ([2], 1.0), + ([1, 2], 1.0), + ([2, 5], 4.0), + ([2, 6], 4.0), + ([5, 6], 4.0), + ([2, 5, 6], 4.0), + ] + def test_expansion_with_blocker(): - st=SimplexTree() - st.insert([0,1],0) - st.insert([0,2],1) - st.insert([0,3],2) - st.insert([1,2],3) - st.insert([1,3],4) - st.insert([2,3],5) - st.insert([2,4],6) - st.insert([3,6],7) - st.insert([4,5],8) - st.insert([4,6],9) - st.insert([5,6],10) - st.insert([6],10) + st = SimplexTree() + st.insert([0, 1], 0) + st.insert([0, 2], 1) + st.insert([0, 3], 2) + st.insert([1, 2], 3) + st.insert([1, 3], 4) + st.insert([2, 3], 5) + st.insert([2, 4], 6) + st.insert([3, 6], 7) + st.insert([4, 5], 8) + st.insert([4, 6], 9) + st.insert([5, 6], 10) + st.insert([6], 10) def blocker(simplex): try: # Block all simplices that contain vertex 6 simplex.index(6) - print(simplex, ' is blocked') + print(simplex, " is blocked") return True except ValueError: - print(simplex, ' is accepted') - st.assign_filtration(simplex, st.filtration(simplex) + 1.) + print(simplex, " is accepted") + st.assign_filtration(simplex, st.filtration(simplex) + 1.0) return False st.expansion_with_blocker(2, blocker) assert st.num_simplices() == 22 assert st.dimension() == 2 - assert st.find([4,5,6]) == False - assert st.filtration([0,1,2]) == 4. - assert st.filtration([0,1,3]) == 5. - assert st.filtration([0,2,3]) == 6. - assert st.filtration([1,2,3]) == 6. + assert st.find([4, 5, 6]) == False + assert st.filtration([0, 1, 2]) == 4.0 + assert st.filtration([0, 1, 3]) == 5.0 + assert st.filtration([0, 2, 3]) == 6.0 + assert st.filtration([1, 2, 3]) == 6.0 st.expansion_with_blocker(3, blocker) assert st.num_simplices() == 23 assert st.dimension() == 3 - assert st.find([4,5,6]) == False - assert st.filtration([0,1,2]) == 4. - assert st.filtration([0,1,3]) == 5. - assert st.filtration([0,2,3]) == 6. - assert st.filtration([1,2,3]) == 6. - assert st.filtration([0,1,2,3]) == 7. + assert st.find([4, 5, 6]) == False + assert st.filtration([0, 1, 2]) == 4.0 + assert st.filtration([0, 1, 3]) == 5.0 + assert st.filtration([0, 2, 3]) == 6.0 + assert st.filtration([1, 2, 3]) == 6.0 + assert st.filtration([0, 1, 2, 3]) == 7.0 diff --git a/src/python/test/test_subsampling.py b/src/python/test/test_subsampling.py index 3431f372..c1cb4e3f 100755 --- a/src/python/test/test_subsampling.py +++ b/src/python/test/test_subsampling.py @@ -16,17 +16,9 @@ __license__ = "MIT" def test_write_off_file_for_tests(): - file = open("subsample.off", "w") - file.write("nOFF\n") - file.write("2 7 0 0\n") - file.write("1.0 1.0\n") - file.write("7.0 0.0\n") - file.write("4.0 6.0\n") - file.write("9.0 6.0\n") - file.write("0.0 14.0\n") - file.write("2.0 19.0\n") - file.write("9.0 17.0\n") - file.close() + gudhi.write_points_to_off_file( + "subsample.off", [[1.0, 1.0], [7.0, 0.0], [4.0, 6.0], [9.0, 6.0], [0.0, 14.0], [2.0, 19.0], [9.0, 17.0]] + ) def test_simple_choose_n_farthest_points_with_a_starting_point(): @@ -34,54 +26,29 @@ def test_simple_choose_n_farthest_points_with_a_starting_point(): i = 0 for point in point_set: # The iteration starts with the given starting point - sub_set = gudhi.choose_n_farthest_points( - points=point_set, nb_points=1, starting_point=i - ) + sub_set = gudhi.choose_n_farthest_points(points=point_set, nb_points=1, starting_point=i) assert sub_set[0] == point_set[i] i = i + 1 # The iteration finds then the farthest - sub_set = gudhi.choose_n_farthest_points( - points=point_set, nb_points=2, starting_point=1 - ) + sub_set = gudhi.choose_n_farthest_points(points=point_set, nb_points=2, starting_point=1) assert sub_set[1] == point_set[3] - sub_set = gudhi.choose_n_farthest_points( - points=point_set, nb_points=2, starting_point=3 - ) + sub_set = gudhi.choose_n_farthest_points(points=point_set, nb_points=2, starting_point=3) assert sub_set[1] == point_set[1] - sub_set = gudhi.choose_n_farthest_points( - points=point_set, nb_points=2, starting_point=0 - ) + sub_set = gudhi.choose_n_farthest_points(points=point_set, nb_points=2, starting_point=0) assert sub_set[1] == point_set[2] - sub_set = gudhi.choose_n_farthest_points( - points=point_set, nb_points=2, starting_point=2 - ) + sub_set = gudhi.choose_n_farthest_points(points=point_set, nb_points=2, starting_point=2) assert sub_set[1] == point_set[0] # Test the limits - assert ( - gudhi.choose_n_farthest_points(points=[], nb_points=0, starting_point=0) == [] - ) - assert ( - gudhi.choose_n_farthest_points(points=[], nb_points=1, starting_point=0) == [] - ) - assert ( - gudhi.choose_n_farthest_points(points=[], nb_points=0, starting_point=1) == [] - ) - assert ( - gudhi.choose_n_farthest_points(points=[], nb_points=1, starting_point=1) == [] - ) + assert gudhi.choose_n_farthest_points(points=[], nb_points=0, starting_point=0) == [] + assert gudhi.choose_n_farthest_points(points=[], nb_points=1, starting_point=0) == [] + assert gudhi.choose_n_farthest_points(points=[], nb_points=0, starting_point=1) == [] + assert gudhi.choose_n_farthest_points(points=[], nb_points=1, starting_point=1) == [] # From off file test for i in range(0, 7): - assert ( - len( - gudhi.choose_n_farthest_points( - off_file="subsample.off", nb_points=i, starting_point=i - ) - ) - == i - ) + assert len(gudhi.choose_n_farthest_points(off_file="subsample.off", nb_points=i, starting_point=i)) == i def test_simple_choose_n_farthest_points_randomed(): @@ -104,10 +71,7 @@ def test_simple_choose_n_farthest_points_randomed(): # From off file test for i in range(0, 7): - assert ( - len(gudhi.choose_n_farthest_points(off_file="subsample.off", nb_points=i)) - == i - ) + assert len(gudhi.choose_n_farthest_points(off_file="subsample.off", nb_points=i)) == i def test_simple_pick_n_random_points(): @@ -130,9 +94,7 @@ def test_simple_pick_n_random_points(): # From off file test for i in range(0, 7): - assert ( - len(gudhi.pick_n_random_points(off_file="subsample.off", nb_points=i)) == i - ) + assert len(gudhi.pick_n_random_points(off_file="subsample.off", nb_points=i)) == i def test_simple_sparsify_points(): @@ -152,31 +114,10 @@ def test_simple_sparsify_points(): ] 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)) - == 7 - ) - assert ( - len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=30.0)) - == 5 - ) - assert ( - 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=89.9)) - == 3 - ) - assert ( - len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=100.0)) - == 2 - ) - assert ( - len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=324.9)) - == 2 - ) - assert ( - len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=325.01)) - == 1 - ) + assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=0.0)) == 7 + assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=30.0)) == 5 + assert 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=89.9)) == 3 + assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=100.0)) == 2 + assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=324.9)) == 2 + assert len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=325.01)) == 1 diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py index 3a004d77..a76b6ce7 100755 --- a/src/python/test/test_wasserstein_distance.py +++ b/src/python/test/test_wasserstein_distance.py @@ -90,10 +90,11 @@ def test_get_essential_parts(): 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) + with pytest.warns(UserWarning): + 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): |