diff options
25 files changed, 369 insertions, 52 deletions
diff --git a/.appveyor.yml b/.appveyor.yml index 5f9c4ef9..125f3cf4 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -26,7 +26,7 @@ environment: PYTHON: "C:\\Python37-x64" - target: Python - CMAKE_FLAGS: -DWITH_GUDHI_EXAMPLE=OFF -DWITH_GUDHI_TEST=OFF -DWITH_GUDHI_UTILITIES=OFF -DWITH_GUDHI_PYTHON=ON -DGMP_INCLUDE_DIR="c:/Tools/vcpkg/installed/x64-windows/include" -DGMP_LIBRARIES="c:/Tools/vcpkg/installed/x64-windows/lib/mpir.lib" -DGMP_LIBRARIES_DIR="c:/Tools/vcpkg/installed/x64-windows/lib" + CMAKE_FLAGS: -DWITH_GUDHI_EXAMPLE=OFF -DWITH_GUDHI_TEST=OFF -DWITH_GUDHI_UTILITIES=OFF -DWITH_GUDHI_PYTHON=ON PYTHON: "C:\\Python37-x64" @@ -48,11 +48,12 @@ install: - pip --version - python -m pip install --upgrade pip - pip install -U setuptools numpy matplotlib scipy Cython pytest + - pip install -U POT build_script: - mkdir build - cd build - - cmake -G "Visual Studio 15 2017 Win64" %CMAKE_FLAGS% -DCMAKE_TOOLCHAIN_FILE=c:/Tools/vcpkg/scripts/buildsystems/vcpkg.cmake .. + - cmake -G "Visual Studio 15 2017 Win64" %CMAKE_FLAGS% -DGMP_INCLUDE_DIR="c:/Tools/vcpkg/installed/x64-windows/include" -DGMP_LIBRARIES="c:/Tools/vcpkg/installed/x64-windows/lib/mpir.lib" -DGMP_LIBRARIES_DIR="c:/Tools/vcpkg/installed/x64-windows/lib" -DCMAKE_TOOLCHAIN_FILE=c:/Tools/vcpkg/scripts/buildsystems/vcpkg.cmake .. - if [%target%]==[Python] ( cd src/python & MSBuild Cython.sln /m /p:Configuration=Release /p:Platform=x64 & diff --git a/.travis.yml b/.travis.yml index bf268057..27514ef5 100644 --- a/.travis.yml +++ b/.travis.yml @@ -56,6 +56,7 @@ before_cache: install: - python3 -m pip install --upgrade pip setuptools wheel - python3 -m pip install --user pytest Cython sphinx sphinxcontrib-bibtex matplotlib numpy scipy + - python3 -m pip install --user POT script: - rm -rf build diff --git a/CMakeGUDHIVersion.txt b/CMakeGUDHIVersion.txt index bc34d9c7..eb2a0666 100644 --- a/CMakeGUDHIVersion.txt +++ b/CMakeGUDHIVersion.txt @@ -1,6 +1,6 @@ set (GUDHI_MAJOR_VERSION 3) set (GUDHI_MINOR_VERSION 0) -set (GUDHI_PATCH_VERSION 0.rc1) +set (GUDHI_PATCH_VERSION 0) set(GUDHI_VERSION ${GUDHI_MAJOR_VERSION}.${GUDHI_MINOR_VERSION}.${GUDHI_PATCH_VERSION}) message(STATUS "GUDHI version : ${GUDHI_VERSION}") diff --git a/Dockerfile_ubuntu b/Dockerfile_for_circleci_image index e149a33a..12f2dc94 100644 --- a/Dockerfile_ubuntu +++ b/Dockerfile_for_circleci_image @@ -51,6 +51,7 @@ RUN pip3 install \ matplotlib \ scipy \ Cython \ + POT \ sphinx \ sphinxcontrib-bibtex diff --git a/Dockerfile_gudhi_installation b/Dockerfile_gudhi_installation new file mode 100644 index 00000000..9fe20730 --- /dev/null +++ b/Dockerfile_gudhi_installation @@ -0,0 +1,65 @@ +FROM ubuntu:19.04 + +# Update and upgrade distribution +RUN apt-get update && \ + apt-get upgrade -y + +# Tools necessary for installing and configuring Ubuntu +RUN apt-get install -y \ + apt-utils \ + locales \ + tzdata + +# Timezone +RUN echo "Europe/Paris" | tee /etc/timezone && \ + ln -fs /usr/share/zoneinfo/Europe/Paris /etc/localtime && \ + dpkg-reconfigure -f noninteractive tzdata + +# Locale with UTF-8 support +RUN echo en_US.UTF-8 UTF-8 >> /etc/locale.gen && \ + locale-gen && \ + update-locale LC_ALL=en_US.UTF-8 LANG=en_US.UTF-8 +ENV LANG en_US.UTF-8 +ENV LANGUAGE en_US:en +ENV LC_ALL en_US.UTF-8 + +# Required for Gudhi compilation +RUN apt-get install -y make \ + g++ \ + cmake \ + graphviz \ + perl \ + texlive-bibtex-extra \ + biber \ + libboost-all-dev \ + libeigen3-dev \ + libgmp3-dev \ + libmpfr-dev \ + libtbb-dev \ + libcgal-dev \ + locales \ + python3 \ + python3-pip \ + python3-pytest \ + python3-tk \ + libfreetype6-dev \ + pkg-config \ + curl + +RUN pip3 install \ + numpy \ + matplotlib \ + scipy \ + Cython + +# apt clean up +RUN apt autoremove && rm -rf /var/lib/apt/lists/* + +RUN curl -LO "https://github.com/GUDHI/gudhi-devel/releases/download/tags%2Fgudhi-release-3.0.0/gudhi.3.0.0.tar.gz" \ +&& tar xf gudhi.3.0.0.tar.gz \ +&& cd gudhi.3.0.0 \ +&& mkdir build && cd build && cmake -DCMAKE_BUILD_TYPE=Release -DWITH_GUDHI_PYTHON=OFF -DPython_ADDITIONAL_VERSIONS=3 .. \ +&& make all test install \ +&& cmake -DWITH_GUDHI_PYTHON=ON . \ +&& cd python \ +&& python3 setup.py install
\ No newline at end of file @@ -2,6 +2,8 @@ [](https://travis-ci.org/GUDHI/gudhi-devel) [](https://circleci.com/gh/GUDHI/gudhi-devel/tree/master) [](https://ci.appveyor.com/project/GUDHI/gudhi-devel/branch/master) +[](https://anaconda.org/conda-forge/gudhi) +[](https://anaconda.org/conda-forge/gudhi)  diff --git a/biblio/bibliography.bib b/biblio/bibliography.bib index d1b2f558..a1b951e0 100644 --- a/biblio/bibliography.bib +++ b/biblio/bibliography.bib @@ -1076,7 +1076,7 @@ language={English} journal = {Computational Geometry: Theory and Applications}, volume = {58}, pages = {70--96}, - doi = "https://doi.org/10.1016/j.comgeo.2016.07.001", + doi = "10.1016/j.comgeo.2016.07.001", year = {2016} } diff --git a/src/Cech_complex/doc/Intro_cech_complex.h b/src/Cech_complex/doc/Intro_cech_complex.h index 90086de7..80c88dc6 100644 --- a/src/Cech_complex/doc/Intro_cech_complex.h +++ b/src/Cech_complex/doc/Intro_cech_complex.h @@ -24,7 +24,7 @@ namespace cech_complex { * \section cechdefinition Čech complex definition * * Čech complex - * <a target="_blank" href="https://en.wikipedia.org/wiki/%C4%8Cech_cohomology">(Wikipedia)</a> is a + * <a target="_blank" href="https://en.wikipedia.org/wiki/%C4%8Cech_complex">(Wikipedia)</a> is a * <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a> constructed * from a proximity graph. The set of all simplices is filtered by the radius of their minimal enclosing ball. * diff --git a/src/Doxyfile.in b/src/Doxyfile.in index 57775498..ec551882 100644 --- a/src/Doxyfile.in +++ b/src/Doxyfile.in @@ -765,7 +765,7 @@ INPUT_ENCODING = UTF-8 # *.md, *.mm, *.dox, *.py, *.f90, *.f, *.for, *.tcl, *.vhd, *.vhdl, *.ucf, # *.qsf, *.as and *.js. -FILE_PATTERNS = +#FILE_PATTERNS = # The RECURSIVE tag can be used to specify whether or not subdirectories should # be searched for input files as well. diff --git a/src/Rips_complex/utilities/ripscomplex.md b/src/Rips_complex/utilities/ripscomplex.md index 03838085..61f31e3c 100644 --- a/src/Rips_complex/utilities/ripscomplex.md +++ b/src/Rips_complex/utilities/ripscomplex.md @@ -99,6 +99,7 @@ where `dim` is the dimension of the homological feature, `birth` and `death` are * `-h [ --help ]` Produce help message * `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. +* `-r [ --max-edge-length ]` (default = inf) Maximal length of an edge for the Rips complex construction. * `-e [ --approximation ]` (default = .5) Epsilon, where the sparse Rips complex is a (1+epsilon)/(1-epsilon)-approximation of the Rips complex. * `-d [ --cpx-dimension ]` (default = INT_MAX) Maximal dimension of the Rips complex we want to compute. * `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. diff --git a/src/Rips_complex/utilities/sparse_rips_persistence.cpp b/src/Rips_complex/utilities/sparse_rips_persistence.cpp index 1a86eafe..cefd8a67 100644 --- a/src/Rips_complex/utilities/sparse_rips_persistence.cpp +++ b/src/Rips_complex/utilities/sparse_rips_persistence.cpp @@ -28,21 +28,24 @@ using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomolog using Point = std::vector<double>; using Points_off_reader = Gudhi::Points_off_reader<Point>; -void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, double& epsilon, +void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, + Filtration_value& threshold, double& epsilon, int& dim_max, int& p, Filtration_value& min_persistence); int main(int argc, char* argv[]) { std::string off_file_points; std::string filediag; + Filtration_value threshold; double epsilon; int dim_max; int p; Filtration_value min_persistence; - program_options(argc, argv, off_file_points, filediag, epsilon, dim_max, p, min_persistence); + program_options(argc, argv, off_file_points, filediag, threshold, epsilon, dim_max, p, min_persistence); Points_off_reader off_reader(off_file_points); - Sparse_rips sparse_rips(off_reader.get_point_cloud(), Gudhi::Euclidean_distance(), epsilon); + Sparse_rips sparse_rips(off_reader.get_point_cloud(), Gudhi::Euclidean_distance(), epsilon, + -std::numeric_limits<Filtration_value>::infinity(), threshold); // Construct the Rips complex in a Simplex Tree Simplex_tree simplex_tree; @@ -73,7 +76,8 @@ int main(int argc, char* argv[]) { return 0; } -void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, double& epsilon, +void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, + Filtration_value& threshold, double& epsilon, int& dim_max, int& p, Filtration_value& min_persistence) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); @@ -84,6 +88,9 @@ void program_options(int argc, char* argv[], std::string& off_file_points, std:: visible.add_options()("help,h", "produce help message")( "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), "Name of file in which the persistence diagram is written. Default print in std::cout")( + "max-edge-length,r", + po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Rips complex construction.")( "approximation,e", po::value<double>(&epsilon)->default_value(.5), "Epsilon, where the sparse Rips complex is a (1+epsilon)-approximation of the Rips complex.")( "cpx-dimension,d", po::value<int>(&dim_max)->default_value(std::numeric_limits<int>::max()), diff --git a/src/cmake/modules/GUDHI_compilation_flags.cmake b/src/cmake/modules/GUDHI_compilation_flags.cmake index 86cd531b..6cd2614d 100644 --- a/src/cmake/modules/GUDHI_compilation_flags.cmake +++ b/src/cmake/modules/GUDHI_compilation_flags.cmake @@ -38,7 +38,7 @@ function(can_cgal_use_cxx11_thread_local) check_cxx_source_compiles("${CGAL_CAN_USE_CXX11_THREAD_LOCAL}" CGAL_CAN_USE_CXX11_THREAD_LOCAL_RESULT) endfunction() -set (CMAKE_CXX_STANDARD 11) +set (CMAKE_CXX_STANDARD 14) enable_testing() diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake index d053e94f..24a34150 100644 --- a/src/cmake/modules/GUDHI_third_party_libraries.cmake +++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake @@ -126,6 +126,7 @@ if( PYTHONINTERP_FOUND ) find_python_module("scipy") find_python_module("sphinx") find_python_module("sklearn") + find_python_module("ot") endif() if(NOT GUDHI_PYTHON_PATH) diff --git a/src/common/doc/installation.h b/src/common/doc/installation.h index 54f86573..2e64bef8 100644 --- a/src/common/doc/installation.h +++ b/src/common/doc/installation.h @@ -5,7 +5,7 @@ * Examples of GUDHI headers inclusion can be found in \ref utilities. * * \section compiling Compiling - * The library uses c++11 and requires <a target="_blank" href="http://www.boost.org/">Boost</a> ≥ 1.56.0 + * The library uses c++14 and requires <a target="_blank" href="http://www.boost.org/">Boost</a> ≥ 1.56.0 * and <a target="_blank" href="https://www.cmake.org/">CMake</a> ≥ 3.1. * It is a multi-platform library and compiles on Linux, Mac OSX and Visual Studio 2015. * diff --git a/src/common/include/gudhi/Unitary_tests_utils.h b/src/common/include/gudhi/Unitary_tests_utils.h index 4ad4dae8..7d039304 100644 --- a/src/common/include/gudhi/Unitary_tests_utils.h +++ b/src/common/include/gudhi/Unitary_tests_utils.h @@ -14,6 +14,7 @@ #include <iostream> #include <limits> // for std::numeric_limits<> +#include <cmath> // for std::fabs template<typename FloatingType > void GUDHI_TEST_FLOAT_EQUALITY_CHECK(FloatingType a, FloatingType b, diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 2756d547..2cc578a6 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -50,6 +50,8 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'euclidean_witness_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'euclidean_strong_witness_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'sktda', ") + # Modules that should not be auto-imported in __init__.py + set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'wasserstein', ") add_gudhi_debug_info("Python version ${PYTHON_VERSION_STRING}") add_gudhi_debug_info("Cython version ${CYTHON_VERSION}") @@ -67,6 +69,8 @@ if(PYTHONINTERP_FOUND) endif() if(SKLEARN_FOUND) add_gudhi_debug_info("Scikit-learn version ${SKLEARN_VERSION}") + if(OT_FOUND) + add_gudhi_debug_info("POT version ${OT_VERSION}") endif() set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DBOOST_RESULT_OF_USE_DECLTYPE', ") @@ -77,7 +81,7 @@ if(PYTHONINTERP_FOUND) if(MSVC) 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++11', ") + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-std=c++14', ") endif(MSVC) if(CMAKE_COMPILER_IS_GNUCXX) set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-frounding-math', ") @@ -204,6 +208,7 @@ if(PYTHONINTERP_FOUND) # Other .py files file(COPY "gudhi/persistence_graphical_tools.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") file(COPY "gudhi/sktda" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/") + file(COPY "gudhi/wasserstein.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi") add_custom_command( OUTPUT gudhi.so @@ -376,37 +381,47 @@ if(PYTHONINTERP_FOUND) # Reader utils add_gudhi_py_test(test_reader_utils) + # Wasserstein + if(OT_FOUND) + add_gudhi_py_test(test_wasserstein_distance) + endif(OT_FOUND) + # Documentation generation is available through sphinx - requires all modules if(SPHINX_PATH) if(MATPLOTLIB_FOUND) if(NUMPY_FOUND) if(SCIPY_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(OT_FOUND) + if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + set (GUDHI_SPHINX_MESSAGE "Generating API documentation with Sphinx in ${CMAKE_CURRENT_BINARY_DIR}/sphinx/") + # User warning - Sphinx is a static pages generator, and configured to work fine with user_version + # Images and biblio warnings because not found on developper version + if (GUDHI_PYTHON_PATH STREQUAL "src/python") + set (GUDHI_SPHINX_MESSAGE "${GUDHI_SPHINX_MESSAGE} \n WARNING : Sphinx is configured for user version, you run it on developper version. Images and biblio will miss") + endif() + # sphinx target requires gudhi.so, because conf.py reads gudhi version from it + add_custom_target(sphinx + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/doc + COMMAND ${CMAKE_COMMAND} -E env "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") + set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") + endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + else(OT_FOUND) + message("++ Python documentation module will not be compiled because POT was not found") set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES") - endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + endif(OT_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") diff --git a/src/python/doc/index.rst b/src/python/doc/index.rst index e379bc23..8f27da0d 100644 --- a/src/python/doc/index.rst +++ b/src/python/doc/index.rst @@ -23,7 +23,7 @@ Alpha complex .. include:: alpha_complex_sum.inc Rips complex -------------- +------------ .. include:: rips_complex_sum.inc @@ -73,6 +73,11 @@ Bottleneck distance .. include:: bottleneck_distance_sum.inc +Wasserstein distance +==================== + +.. include:: wasserstein_distance_sum.inc + Persistence graphical tools =========================== diff --git a/src/python/doc/installation.rst b/src/python/doc/installation.rst index d8b6f861..7699a5bb 100644 --- a/src/python/doc/installation.rst +++ b/src/python/doc/installation.rst @@ -8,11 +8,11 @@ Installation Conda ***** The easiest way to install the Python version of GUDHI is using -`conda <https://gudhi.inria.fr/licensing/>`_. +`conda <https://gudhi.inria.fr/conda/>`_. Compiling ********* -The library uses c++11 and requires `Boost <https://www.boost.org/>`_ ≥ 1.56.0, +The library uses c++14 and requires `Boost <https://www.boost.org/>`_ ≥ 1.56.0, `CMake <https://www.cmake.org/>`_ ≥ 3.1 to generate makefiles, `NumPy <http://numpy.org>`_ and `Cython <https://www.cython.org/>`_ to compile the GUDHI Python module. @@ -138,7 +138,7 @@ Documentation To build the documentation, `sphinx-doc <http://www.sphinx-doc.org>`_ and `sphinxcontrib-bibtex <https://sphinxcontrib-bibtex.readthedocs.io>`_ are -required. As the documentation is auto-tested, `CGAL`_, `Eigen3`_, +required. As the documentation is auto-tested, `CGAL`_, `Eigen`_, `Matplotlib`_, `NumPy`_ and `SciPy`_ are also mandatory to build the documentation. @@ -215,12 +215,20 @@ 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>` +Python Optimal Transport +======================== + +The :doc:`Wasserstein distance </wasserstein_distance_user>` +module requires `POT <https://pot.readthedocs.io/>`_, a library that provides +several solvers for optimization problems related to Optimal Transport. + SciPy ===== -The :doc:`persistence graphical tools </persistence_graphical_tools_user>` -module requires `SciPy <http://scipy.org>`_, a Python-based ecosystem of -open-source software for mathematics, science, and engineering. +The :doc:`persistence graphical tools </persistence_graphical_tools_user>` and +:doc:`Wasserstein distance </wasserstein_distance_user>` modules require `SciPy +<http://scipy.org>`_, a Python-based ecosystem of open-source software for +mathematics, science, and engineering. Threading Building Blocks ========================= diff --git a/src/python/doc/wasserstein_distance_sum.inc b/src/python/doc/wasserstein_distance_sum.inc new file mode 100644 index 00000000..ffd4d312 --- /dev/null +++ b/src/python/doc/wasserstein_distance_sum.inc @@ -0,0 +1,14 @@ +.. table:: + :widths: 30 50 20 + + +-----------------------------------------------------------------+----------------------------------------------------------------------+------------------------------------------------------------------+ + | .. figure:: | The p-Wasserstein distance measures the similarity between two | :Author: Theo Lacombe | + | ../../doc/Bottleneck_distance/perturb_pd.png | persistence diagrams. It's the minimum value c that can be achieved | | + | :figclass: align-center | by a perfect matching between the points of the two diagrams (+ all | :Introduced in: GUDHI 3.1.0 | + | | diagonal points), where the value of a matching is defined as the | | + | Wasserstein distance is the p-th root of the sum of the | p-th root of the sum of all edge lengths to the power p. Edge lengths| :Copyright: MIT | + | edge lengths to the power p. | are measured in norm q, for :math:`1 \leq q \leq \infty`. | | + | | | :Requires: Python Optimal Transport (POT) :math:`\geq` 0.5.1 | + +-----------------------------------------------------------------+----------------------------------------------------------------------+------------------------------------------------------------------+ + | * :doc:`wasserstein_distance_user` | | + +-----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/doc/wasserstein_distance_user.rst b/src/python/doc/wasserstein_distance_user.rst new file mode 100644 index 00000000..a049cfb5 --- /dev/null +++ b/src/python/doc/wasserstein_distance_user.rst @@ -0,0 +1,40 @@ +:orphan: + +.. To get rid of WARNING: document isn't included in any toctree + +Wasserstein distance user manual +================================ +Definition +---------- + +.. include:: wasserstein_distance_sum.inc + +This implementation is based on ideas from "Large Scale Computation of Means and Cluster for Persistence Diagrams via Optimal Transport". + +Function +-------- +.. autofunction:: gudhi.wasserstein.wasserstein_distance + + +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. + +.. testcode:: + + import gudhi.wasserstein + import numpy as np + + diag1 = np.array([[2.7, 3.7],[9.6, 14.],[34.2, 34.974]]) + diag2 = np.array([[2.8, 4.45],[9.5, 14.1]]) + + message = "Wasserstein distance value = " + '%.2f' % gudhi.wasserstein.wasserstein_distance(diag1, diag2, q=2., p=1.) + print(message) + +The output is: + +.. testoutput:: + + Wasserstein distance value = 1.45 diff --git a/src/python/gudhi/__init__.py.in b/src/python/gudhi/__init__.py.in index 28bab0e1..02888fff 100644 --- a/src/python/gudhi/__init__.py.in +++ b/src/python/gudhi/__init__.py.in @@ -21,13 +21,16 @@ __debug_info__ = @GUDHI_PYTHON_DEBUG_INFO@ from sys import exc_info from importlib import import_module -__all__ = [@GUDHI_PYTHON_MODULES@] +__all__ = [@GUDHI_PYTHON_MODULES@ @GUDHI_PYTHON_MODULES_EXTRA@] __available_modules = '' __missing_modules = '' -# try to import * from gudhi.__module_name -for __module_name in __all__: +# 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 +# speed up the basic import) +for __module_name in [@GUDHI_PYTHON_MODULES@]: try: __module = import_module('gudhi.' + __module_name) try: diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 9f490271..90ddf6dd 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -362,7 +362,7 @@ cdef class SimplexTree: value than its faces by increasing the filtration values. :returns: True if any filtration value was modified, - False if the filtration was already non-decreasing. + False if the filtration was already non-decreasing. :rtype: bool diff --git a/src/python/gudhi/wasserstein.py b/src/python/gudhi/wasserstein.py new file mode 100644 index 00000000..eba7c6d5 --- /dev/null +++ b/src/python/gudhi/wasserstein.py @@ -0,0 +1,99 @@ +import numpy as np +import scipy.spatial.distance as sc +try: + import ot +except ImportError: + print("POT (Python Optimal Transport) package is not installed. Try to run $ conda install -c conda-forge pot ; or $ pip install POT") + +""" 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): Theo Lacombe + + Copyright (C) 2019 Inria + + Modification(s): + - YYYY/MM Author: Description of the modification +""" + +def _proj_on_diag(X): + ''' + :param X: (n x 2) array encoding the points of a persistent diagram. + :returns: (n x 2) array encoding the (respective orthogonal) projections of the points onto the diagonal + ''' + Z = (X[:,0] + X[:,1]) / 2. + return np.array([Z , Z]).T + + +def _build_dist_matrix(X, Y, p=2., q=2.): + ''' + :param X: (n x 2) numpy.array encoding the (points of the) first diagram. + :param Y: (m x 2) numpy.array encoding the second diagram. + :param q: Ground metric (i.e. norm l_q). + :param p: exponent for the Wasserstein metric. + :returns: (n+1) x (m+1) np.array encoding the cost matrix C. + For 1 <= i <= n, 1 <= j <= m, C[i,j] encodes the distance between X[i] and Y[j], while C[i, m+1] (resp. C[n+1, j]) encodes the distance (to the p) between X[i] (resp Y[j]) and its orthogonal proj onto the diagonal. + note also that C[n+1, m+1] = 0 (it costs nothing to move from the diagonal to the diagonal). + ''' + Xdiag = _proj_on_diag(X) + Ydiag = _proj_on_diag(Y) + if np.isinf(q): + C = sc.cdist(X,Y, metric='chebyshev')**p + Cxd = np.linalg.norm(X - Xdiag, ord=q, axis=1)**p + Cdy = np.linalg.norm(Y - Ydiag, ord=q, axis=1)**p + else: + C = sc.cdist(X,Y, metric='minkowski', p=q)**p + Cxd = np.linalg.norm(X - Xdiag, ord=q, axis=1)**p + Cdy = np.linalg.norm(Y - Ydiag, ord=q, axis=1)**p + Cf = np.hstack((C, Cxd[:,None])) + Cdy = np.append(Cdy, 0) + + Cf = np.vstack((Cf, Cdy[None,:])) + + return Cf + + +def _perstot(X, p, q): + ''' + :param X: (n x 2) numpy.array (points of a given diagram). + :param q: Ground metric on the (upper-half) plane (i.e. norm l_q in R^2); Default value is 2 (Euclidean norm). + :param p: exponent for Wasserstein; Default value is 2. + :returns: float, the total persistence of the diagram (that is, its distance to the empty diagram). + ''' + Xdiag = _proj_on_diag(X) + return (np.sum(np.linalg.norm(X - Xdiag, ord=q, axis=1)**p))**(1./p) + + +def wasserstein_distance(X, Y, p=2., q=2.): + ''' + :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 q: Ground metric on the (upper-half) plane (i.e. norm l_q in R^2); Default value is 2 (euclidean norm). + :param p: exponent for Wasserstein; Default value is 2. + :returns: the p-Wasserstein distance (1 <= p < infinity) with respect to the q-norm as ground metric. + :rtype: float + ''' + n = len(X) + m = len(Y) + + # handle empty diagrams + if X.size == 0: + if Y.size == 0: + return 0. + else: + return _perstot(Y, p, q) + elif Y.size == 0: + return _perstot(X, p, q) + + M = _build_dist_matrix(X, Y, p=p, q=q) + a = np.full(n+1, 1. / (n + m) ) # weight vector of the input diagram. Uniform here. + a[-1] = a[-1] * m # normalized so that we have a probability measure, required by POT + b = np.full(m+1, 1. / (n + m) ) # weight vector of the input diagram. Uniform here. + b[-1] = b[-1] * n # so that we have a probability measure, required by POT + + # Comptuation of the otcost using the ot.emd2 library. + # Note: it is the squared Wasserstein distance. + # The default numItermax=100000 is not sufficient for some examples with 5000 points, what is a good value? + ot_cost = (n+m) * ot.emd2(a, b, M, numItermax=2000000) + + return ot_cost ** (1./p) + diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index b3553d32..96353cc4 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -15,6 +15,8 @@ #include <gudhi/Alpha_complex.h> #include <CGAL/Epick_d.h> +#include <boost/range/adaptor/transformed.hpp> + #include "Simplex_tree_interface.h" #include <iostream> @@ -31,7 +33,10 @@ class Alpha_complex_interface { public: Alpha_complex_interface(const std::vector<std::vector<double>>& points) { - alpha_complex_ = new Alpha_complex<Dynamic_kernel>(points); + auto mkpt = [](std::vector<double> const& vec){ + return Point_d(vec.size(), vec.begin(), vec.end()); + }; + alpha_complex_ = new Alpha_complex<Dynamic_kernel>(boost::adaptors::transform(points, mkpt)); } Alpha_complex_interface(const std::string& off_file_name, bool from_file = true) { @@ -45,9 +50,9 @@ class Alpha_complex_interface { std::vector<double> get_point(int vh) { std::vector<double> vd; try { - Point_d ph = alpha_complex_->get_point(vh); + Point_d const& ph = alpha_complex_->get_point(vh); for (auto coord = ph.cartesian_begin(); coord < ph.cartesian_end(); coord++) - vd.push_back(*coord); + vd.push_back(CGAL::to_double(*coord)); } catch (std::out_of_range const&) { // std::out_of_range is thrown in case not found. Other exceptions must be re-thrown } diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py new file mode 100755 index 00000000..a6bf9901 --- /dev/null +++ b/src/python/test/test_wasserstein_distance.py @@ -0,0 +1,48 @@ +from gudhi.wasserstein import wasserstein_distance +import numpy as np + +""" 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): Theo Lacombe + + Copyright (C) 2019 Inria + + Modification(s): + - YYYY/MM Author: Description of the modification +""" + +__author__ = "Theo Lacombe" +__copyright__ = "Copyright (C) 2019 Inria" +__license__ = "MIT" + + +def test_basic_wasserstein(): + 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]]) + diag3 = np.array([[0, 2], [4, 6]]) + diag4 = np.array([[0, 3], [4, 8]]) + emptydiag = np.array([[]]) + + assert wasserstein_distance(emptydiag, emptydiag, q=2., p=1.) == 0. + assert wasserstein_distance(emptydiag, emptydiag, q=np.inf, p=1.) == 0. + assert wasserstein_distance(emptydiag, emptydiag, q=np.inf, p=2.) == 0. + assert wasserstein_distance(emptydiag, emptydiag, q=2., p=2.) == 0. + + assert wasserstein_distance(diag3, emptydiag, q=np.inf, p=1.) == 2. + assert wasserstein_distance(diag3, emptydiag, q=1., p=1.) == 4. + + assert wasserstein_distance(diag4, emptydiag, q=1., p=2.) == 5. # thank you Pythagorician triplets + assert wasserstein_distance(diag4, emptydiag, q=np.inf, p=2.) == 2.5 + assert wasserstein_distance(diag4, emptydiag, q=2., p=2.) == 3.5355339059327378 + + assert wasserstein_distance(diag1, diag2, q=2., p=1.) == 1.4453593023967701 + assert wasserstein_distance(diag1, diag2, q=2.35, p=1.74) == 0.9772734057168739 + + assert wasserstein_distance(diag1, emptydiag, q=2.35, p=1.7863) == 3.141592214572228 + + assert wasserstein_distance(diag3, diag4, q=1., p=1.) == 3. + assert wasserstein_distance(diag3, diag4, q=np.inf, p=1.) == 3. # no diag matching here + assert wasserstein_distance(diag3, diag4, q=np.inf, p=2.) == np.sqrt(5) + assert wasserstein_distance(diag3, diag4, q=1., p=2.) == np.sqrt(5) + assert wasserstein_distance(diag3, diag4, q=4.5, p=2.) == np.sqrt(5) + |