diff options
Diffstat (limited to 'src/python')
-rw-r--r-- | src/python/CMakeLists.txt | 26 | ||||
-rw-r--r-- | src/python/doc/alpha_complex_user.rst | 7 | ||||
-rw-r--r-- | src/python/doc/installation.rst | 99 | ||||
-rw-r--r-- | src/python/doc/nerve_gic_complex_user.rst | 2 | ||||
-rw-r--r-- | src/python/doc/persistence_graphical_tools_user.rst | 2 | ||||
-rw-r--r-- | src/python/gudhi/__init__.py.in | 4 | ||||
-rw-r--r-- | src/python/gudhi/alpha_complex.pyx | 29 | ||||
-rw-r--r-- | src/python/gudhi/persistence_graphical_tools.py | 349 | ||||
-rw-r--r-- | src/python/gudhi/simplex_tree.pxd | 5 | ||||
-rw-r--r-- | src/python/gudhi/simplex_tree.pyx | 86 | ||||
-rw-r--r-- | src/python/gudhi/weighted_rips_complex.py | 6 | ||||
-rw-r--r-- | src/python/include/Alpha_complex_interface.h | 10 | ||||
-rw-r--r-- | src/python/include/Simplex_tree_interface.h | 16 | ||||
-rw-r--r-- | src/python/pyproject.toml | 2 | ||||
-rw-r--r-- | src/python/setup.py.in | 4 | ||||
-rwxr-xr-x | src/python/test/test_alpha_complex.py | 27 | ||||
-rw-r--r-- | src/python/test/test_persistence_graphical_tools.py | 121 | ||||
-rwxr-xr-x | src/python/test/test_simplex_tree.py | 131 |
18 files changed, 652 insertions, 274 deletions
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index ab5bc56f..4b04c487 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -148,10 +148,6 @@ if(PYTHONINTERP_FOUND) add_gudhi_debug_info("Eigen3 version ${EIGEN3_VERSION}") # No problem, even if no CGAL found set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_EIGEN3_ENABLED', ") - set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DGUDHI_USE_EIGEN3', ") - set(GUDHI_USE_EIGEN3 "True") - else (EIGEN3_FOUND) - set(GUDHI_USE_EIGEN3 "False") endif (EIGEN3_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_reader', ") @@ -180,6 +176,16 @@ if(PYTHONINTERP_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") endif () + # from windows vcpkg eigen 3.4.0#2 : build fails with + # error C2440: '<function-style-cast>': cannot convert from 'Eigen::EigenBase<Derived>::Index' to '__gmp_expr<mpq_t,mpq_t>' + # cf. https://gitlab.com/libeigen/eigen/-/issues/2476 + # Workaround is to compile with '-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int' + if (FORCE_EIGEN_DEFAULT_DENSE_INDEX_TYPE_TO_INT) + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int', ") + endif() + + + add_gudhi_debug_info("Boost version ${Boost_VERSION}") if(CGAL_FOUND) if(NOT CGAL_VERSION VERSION_LESS 5.3.0) # CGAL_HEADER_ONLY has been dropped for CGAL >= 5.3. Only the header-only version is supported. @@ -214,13 +220,14 @@ if(PYTHONINTERP_FOUND) endif(NOT GMP_LIBRARIES_DIR) add_GUDHI_PYTHON_lib_dir(${GMP_LIBRARIES_DIR}) message("** Add gmp ${GMP_LIBRARIES_DIR}") + # When FORCE_CGAL_NOT_TO_BUILD_WITH_GMPXX is set, not defining CGAL_USE_GMPXX is sufficient enough if(GMPXX_FOUND) add_gudhi_debug_info("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}") set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_GMPXX', ") add_GUDHI_PYTHON_lib("${GMPXX_LIBRARIES}") add_GUDHI_PYTHON_lib_dir(${GMPXX_LIBRARIES_DIR}) message("** Add gmpxx ${GMPXX_LIBRARIES_DIR}") - endif(GMPXX_FOUND) + endif() endif(GMP_FOUND) if(MPFR_FOUND) add_gudhi_debug_info("MPFR_LIBRARIES = ${MPFR_LIBRARIES}") @@ -264,6 +271,10 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_INCLUDE_DIRS "${GUDHI_PYTHON_INCLUDE_DIRS}'${TBB_INCLUDE_DIRS}', ") endif() + if(DEBUG_TRACES) + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DDEBUG_TRACES', ") + endif(DEBUG_TRACES) + if(UNIX AND WITH_GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS) set( GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}") endif(UNIX AND WITH_GUDHI_PYTHON_RUNTIME_LIBRARY_DIRS) @@ -591,6 +602,11 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_sklearn_cubical_persistence) endif() + # persistence graphical tools + if(MATPLOTLIB_FOUND) + add_gudhi_py_test(test_persistence_graphical_tools) + endif() + # Set missing or not modules set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES") else(PYBIND11_FOUND AND CYTHON_FOUND) diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst index cfd22742..9e67d38a 100644 --- a/src/python/doc/alpha_complex_user.rst +++ b/src/python/doc/alpha_complex_user.rst @@ -27,7 +27,8 @@ Remarks If you pass :code:`precision = 'exact'` to :func:`~gudhi.AlphaComplex.__init__`, the filtration values are the exact ones converted to float. This can be very slow. If you pass :code:`precision = 'safe'` (the default), the filtration values are only - guaranteed to have a small multiplicative error compared to the exact value. + guaranteed to have a small multiplicative error compared to the exact value, see + :func:`~gudhi.AlphaComplex.set_float_relative_precision` to modify the precision. A drawback, when computing persistence, is that an empty exact interval [10^12,10^12] may become a non-empty approximate interval [10^12,10^12+10^6]. Using :code:`precision = 'fast'` makes the computations slightly faster, and the combinatorics are still exact, but @@ -177,11 +178,11 @@ Weighted version ^^^^^^^^^^^^^^^^ A weighted version for Alpha complex is available. It is like a usual Alpha complex, but based on a -`CGAL regular triangulation <https://doc.cgal.org/latest/Triangulation/index.html#title20>`_. +`CGAL regular triangulation <https://doc.cgal.org/latest/Triangulation/index.html#TriangulationSecRT>`_. This example builds the weighted alpha-complex of a small molecule, where atoms have different sizes. It is taken from -`CGAL 3d weighted alpha shapes <https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title13>`_. +`CGAL 3d weighted alpha shapes <https://doc.cgal.org/latest/Alpha_shapes_3/index.html#AlphaShape_3DExampleforWeightedAlphaShapes>`_. Then, it is asked to display information about the alpha complex. diff --git a/src/python/doc/installation.rst b/src/python/doc/installation.rst index 35c344e3..cff84691 100644 --- a/src/python/doc/installation.rst +++ b/src/python/doc/installation.rst @@ -33,25 +33,19 @@ Compiling These instructions are for people who want to compile gudhi from source, they are unnecessary if you installed a binary package of Gudhi as above. They assume that you have downloaded a `release <https://github.com/GUDHI/gudhi-devel/releases>`_, -with a name like `gudhi.3.2.0.tar.gz`, then run `tar xf gudhi.3.2.0.tar.gz`, which -created a directory `gudhi.3.2.0`, hereinafter referred to as `/path-to-gudhi/`. +with a name like `gudhi.3.X.Y.tar.gz`, then run `tar xf gudhi.3.X.Y.tar.gz`, which +created a directory `gudhi.3.X.Y`, hereinafter referred to as `/path-to-gudhi/`. If you are instead using a git checkout, beware that the paths are a bit different, and in particular the `python/` subdirectory is actually `src/python/` there. -The library uses c++14 and requires `Boost <https://www.boost.org/>`_ :math:`\geq` 1.56.0, +The library uses c++14 and requires `Boost <https://www.boost.org/>`_ :math:`\geq` 1.66.0, `CMake <https://www.cmake.org/>`_ :math:`\geq` 3.5 to generate makefiles, -`NumPy <http://numpy.org>`_ :math:`\geq` 1.15.0, `Cython <https://www.cython.org/>`_ and -`pybind11 <https://github.com/pybind/pybind11>`_ to compile -the GUDHI Python module. -It is a multi-platform library and compiles on Linux, Mac OSX and Visual -Studio 2017 or later. +Python :math:`\geq` 3.5, `NumPy <http://numpy.org>`_ :math:`\geq` 1.15.0, `Cython <https://www.cython.org/>`_ +:math:`\geq` 0.27 and `pybind11 <https://github.com/pybind/pybind11>`_ to compile the GUDHI Python module. +It is a multi-platform library and compiles on Linux, Mac OSX and Visual Studio 2017 or later. -On `Windows <https://wiki.python.org/moin/WindowsCompilers>`_ , only Python -:math:`\geq` 3.5 are available because of the required Visual Studio version. - -On other systems, if you have several Python/python installed, the version 2.X -will be used by default, but you can force it by adding +If you have several Python/python installed, the version 2.X may be used by default, but you can force it by adding :code:`-DPython_ADDITIONAL_VERSIONS=3` to the cmake command. GUDHI Python module compilation @@ -142,54 +136,63 @@ If :code:`import gudhi` succeeds, please have a look to debug information: .. code-block:: python - import gudhi - print(gudhi.__debug_info__) + import gudhi as gd + print(gd.__debug_info__) + print("+ Installed modules are: " + gd.__available_modules) + print("+ Missing modules are: " + gd.__missing_modules) You shall have something like: .. code-block:: none - Python version 2.7.15 - Cython version 0.26.1 - Numpy version 1.14.1 - Eigen3 version 3.1.1 - Installed modules are: off_reader;simplex_tree;rips_complex; - cubical_complex;periodic_cubical_complex;reader_utils;witness_complex; - strong_witness_complex;alpha_complex; - Missing modules are: bottleneck_distance;nerve_gic;subsampling; - tangential_complex;persistence_graphical_tools; - euclidean_witness_complex;euclidean_strong_witness_complex; - CGAL version 4.7.1000 - GMP_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmp.so - GMPXX_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmpxx.so - TBB version 9107 found and used + Pybind11 version 2.8.1 + Python version 3.7.12 + Cython version 0.29.25 + Numpy version 1.21.4 + Boost version 1.77.0 + + Installed modules are: off_reader;simplex_tree;rips_complex;cubical_complex;periodic_cubical_complex; + persistence_graphical_tools;reader_utils;witness_complex;strong_witness_complex; + + Missing modules are: bottleneck;nerve_gic;subsampling;tangential_complex;alpha_complex;euclidean_witness_complex; + euclidean_strong_witness_complex; -Here, you can see that bottleneck_distance, nerve_gic, subsampling and -tangential_complex are missing because of the CGAL version. -persistence_graphical_tools is not available as matplotlib is not -available. +Here, you can see that the modules that need CGAL are missing, because CGAL is not installed. +:code:`persistence_graphical_tools` is installed, but +`its functions <https://gudhi.inria.fr/python/latest/persistence_graphical_tools_ref.html>`_ will produce an error as +matplotlib is not available. Unitary tests cannot be run as pytest is missing. A complete configuration would be : .. code-block:: none - Python version 3.6.5 - Cython version 0.28.2 - Pytest version 3.3.2 - Matplotlib version 2.2.2 - Numpy version 1.14.5 - Eigen3 version 3.3.4 - Installed modules are: off_reader;simplex_tree;rips_complex; - cubical_complex;periodic_cubical_complex;persistence_graphical_tools; - reader_utils;witness_complex;strong_witness_complex; - persistence_graphical_tools;bottleneck_distance;nerve_gic;subsampling; - tangential_complex;alpha_complex;euclidean_witness_complex; - euclidean_strong_witness_complex; - CGAL header only version 4.11.0 + Pybind11 version 2.8.1 + Python version 3.9.7 + Cython version 0.29.24 + Pytest version 6.2.5 + Matplotlib version 3.5.0 + Numpy version 1.21.4 + Scipy version 1.7.3 + Scikit-learn version 1.0.1 + POT version 0.8.0 + HNSWlib found + PyKeOps version [pyKeOps]: 1.5 + EagerPy version 0.30.0 + TensorFlow version 2.7.0 + Sphinx version 4.3.0 + Sphinx-paramlinks version 0.5.2 + python_docs_theme found + Eigen3 version 3.4.0 + Boost version 1.74.0 + CGAL version 5.3 GMP_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmp.so GMPXX_LIBRARIES = /usr/lib/x86_64-linux-gnu/libgmpxx.so + MPFR_LIBRARIES = /usr/lib/x86_64-linux-gnu/libmpfr.so TBB version 9107 found and used + + Installed modules are: bottleneck;off_reader;simplex_tree;rips_complex;cubical_complex;periodic_cubical_complex; + persistence_graphical_tools;reader_utils;witness_complex;strong_witness_complex;nerve_gic;subsampling; + tangential_complex;alpha_complex;euclidean_witness_complex;euclidean_strong_witness_complex; + + Missing modules are: + Documentation ============= @@ -345,8 +348,8 @@ You can still deactivate LaTeX rendering by saying: .. code-block:: python - import gudhi - gudhi.persistence_graphical_tools._gudhi_matplotlib_use_tex=False + import gudhi as gd + gd.persistence_graphical_tools._gudhi_matplotlib_use_tex=False PyKeOps ------- diff --git a/src/python/doc/nerve_gic_complex_user.rst b/src/python/doc/nerve_gic_complex_user.rst index 0b820abf..8633cadb 100644 --- a/src/python/doc/nerve_gic_complex_user.rst +++ b/src/python/doc/nerve_gic_complex_user.rst @@ -12,7 +12,7 @@ Definition Visualizations of the simplicial complexes can be done with either neato (from `graphviz <http://www.graphviz.org/>`_), `geomview <http://www.geomview.org/>`_, -`KeplerMapper <https://github.com/MLWave/kepler-mapper>`_. +`KeplerMapper <https://github.com/scikit-tda/kepler-mapper>`_. Input point clouds are assumed to be OFF files (cf. `OFF file format <fileformats.html#off-file-format>`_). Covers diff --git a/src/python/doc/persistence_graphical_tools_user.rst b/src/python/doc/persistence_graphical_tools_user.rst index d95b9d2b..e1d28c71 100644 --- a/src/python/doc/persistence_graphical_tools_user.rst +++ b/src/python/doc/persistence_graphical_tools_user.rst @@ -60,7 +60,7 @@ of shape (N x 2) encoding a persistence diagram (in a given dimension). import matplotlib.pyplot as plt import gudhi import numpy as np - d = np.array([[0, 1], [1, 2], [1, np.inf]]) + d = np.array([[0., 1.], [1., 2.], [1., np.inf]]) gudhi.plot_persistence_diagram(d) plt.show() diff --git a/src/python/gudhi/__init__.py.in b/src/python/gudhi/__init__.py.in index 3043201a..79e12fbc 100644 --- a/src/python/gudhi/__init__.py.in +++ b/src/python/gudhi/__init__.py.in @@ -23,10 +23,6 @@ __all__ = [@GUDHI_PYTHON_MODULES@ @GUDHI_PYTHON_MODULES_EXTRA@] __available_modules = '' __missing_modules = '' -# For unitary tests purpose -# could use "if 'collapse_edges' in gudhi.__all__" when collapse edges will have a python module -__GUDHI_USE_EIGEN3 = @GUDHI_USE_EIGEN3@ - # Try to import * from gudhi.__module_name for default modules. # Extra modules require an explicit import by the user (mostly because of # unusual dependencies, but also to avoid cluttering namespace gudhi and diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx index a4888914..375e1561 100644 --- a/src/python/gudhi/alpha_complex.pyx +++ b/src/python/gudhi/alpha_complex.pyx @@ -31,6 +31,10 @@ cdef extern from "Alpha_complex_interface.h" namespace "Gudhi": Alpha_complex_interface(vector[vector[double]] points, vector[double] weights, bool fast_version, bool exact_version) nogil except + vector[double] get_point(int vertex) nogil except + void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool default_filtration_value) nogil except + + @staticmethod + void set_float_relative_precision(double precision) nogil + @staticmethod + double get_float_relative_precision() nogil # AlphaComplex python interface cdef class AlphaComplex: @@ -133,3 +137,28 @@ cdef class AlphaComplex: self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr, mas, compute_filtration) return stree + + @staticmethod + def set_float_relative_precision(precision): + """ + :param precision: When the AlphaComplex is constructed with :code:`precision = 'safe'` (the default), + one can set the float relative precision of filtration values computed in + :func:`~gudhi.AlphaComplex.create_simplex_tree`. + Default is :code:`1e-5` (cf. :func:`~gudhi.AlphaComplex.get_float_relative_precision`). + For more details, please refer to + `CGAL::Lazy_exact_nt<NT>::set_relative_precision_of_to_double <https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html>`_ + :type precision: float + """ + if precision <=0. or precision >= 1.: + raise ValueError("Relative precision value must be strictly greater than 0 and strictly lower than 1") + Alpha_complex_interface.set_float_relative_precision(precision) + + @staticmethod + def get_float_relative_precision(): + """ + :returns: The float relative precision of filtration values computation in + :func:`~gudhi.AlphaComplex.create_simplex_tree` when the AlphaComplex is constructed with + :code:`precision = 'safe'` (the default). + :rtype: float + """ + return Alpha_complex_interface.get_float_relative_precision() diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py index 848dc03e..7ed11360 100644 --- a/src/python/gudhi/persistence_graphical_tools.py +++ b/src/python/gudhi/persistence_graphical_tools.py @@ -12,6 +12,9 @@ from os import path from math import isfinite import numpy as np from functools import lru_cache +import warnings +import errno +import os from gudhi.reader_utils import read_persistence_intervals_in_dimension from gudhi.reader_utils import read_persistence_intervals_grouped_by_dimension @@ -22,6 +25,7 @@ __license__ = "MIT" _gudhi_matplotlib_use_tex = True + def __min_birth_max_death(persistence, band=0.0): """This function returns (min_birth, max_death) from the persistence. @@ -44,20 +48,46 @@ def __min_birth_max_death(persistence, band=0.0): min_birth = float(interval[1][0]) if band > 0.0: max_death += band + # can happen if only points at inf death + if min_birth == max_death: + max_death = max_death + 1.0 return (min_birth, max_death) def _array_handler(a): - ''' + """ :param a: if array, assumes it is a (n x 2) np.array and return a persistence-compatible list (padding with 0), so that the plot can be performed seamlessly. - ''' - if isinstance(a[0][1], np.float64) or isinstance(a[0][1], float): + """ + if isinstance(a[0][1], (np.floating, float)): return [[0, x] for x in a] else: return a + +def _limit_to_max_intervals(persistence, max_intervals, key): + """This function returns truncated persistence if length is bigger than max_intervals. + :param persistence: Persistence intervals values list. Can be grouped by dimension or not. + :type persistence: an array of (dimension, array of (birth, death)) or an array of (birth, death). + :param max_intervals: maximal number of intervals to display. + Selected intervals are those with the longest life time. Set it + to 0 to see all. Default value is 1000. + :type max_intervals: int. + :param key: key function for sort algorithm. + :type key: function or lambda. + """ + if max_intervals > 0 and max_intervals < len(persistence): + warnings.warn( + "There are %s intervals given as input, whereas max_intervals is set to %s." + % (len(persistence), max_intervals) + ) + # Sort by life time, then takes only the max_intervals elements + return sorted(persistence, key=key, reverse=True)[:max_intervals] + else: + return persistence + + @lru_cache(maxsize=1) def _matplotlib_can_use_tex(): """This function returns True if matplotlib can deal with LaTeX, False otherwise. @@ -65,17 +95,17 @@ def _matplotlib_can_use_tex(): """ try: from matplotlib import checkdep_usetex + return checkdep_usetex(True) - except ImportError: - print("This function is not available, you may be missing matplotlib.") + except ImportError as import_error: + warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.") def plot_persistence_barcode( persistence=[], persistence_file="", alpha=0.6, - max_intervals=1000, - max_barcodes=1000, + max_intervals=20000, inf_delta=0.1, legend=False, colormap=None, @@ -97,7 +127,7 @@ def plot_persistence_barcode( :type alpha: float. :param max_intervals: maximal number of intervals to display. Selected intervals are those with the longest life time. Set it - to 0 to see all. Default value is 1000. + to 0 to see all. Default value is 20000. :type max_intervals: int. :param inf_delta: Infinity is placed at :code:`((max_death - min_birth) x inf_delta)` above :code:`max_death` value. A reasonable value is @@ -119,99 +149,68 @@ def plot_persistence_barcode( import matplotlib.pyplot as plt import matplotlib.patches as mpatches from matplotlib import rc + if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex(): - plt.rc('text', usetex=True) - plt.rc('font', family='serif') + plt.rc("text", usetex=True) + plt.rc("font", family="serif") else: - plt.rc('text', usetex=False) - plt.rc('font', family='DejaVu Sans') + plt.rc("text", usetex=False) + plt.rc("font", family="DejaVu Sans") if persistence_file != "": if path.isfile(persistence_file): # Reset persistence persistence = [] - diag = read_persistence_intervals_grouped_by_dimension( - persistence_file=persistence_file - ) + diag = read_persistence_intervals_grouped_by_dimension(persistence_file=persistence_file) for key in diag.keys(): for persistence_interval in diag[key]: persistence.append((key, persistence_interval)) else: - print("file " + persistence_file + " not found.") - return None - - persistence = _array_handler(persistence) - - if max_barcodes != 1000: - print("Deprecated parameter. It has been replaced by max_intervals") - max_intervals = max_barcodes - - if max_intervals > 0 and max_intervals < len(persistence): - # Sort by life time, then takes only the max_intervals elements - persistence = sorted( - persistence, - key=lambda life_time: life_time[1][1] - life_time[1][0], - reverse=True, - )[:max_intervals] - - if colormap == None: - colormap = plt.cm.Set1.colors - if axes == None: - fig, axes = plt.subplots(1, 1) + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file) - persistence = sorted(persistence, key=lambda birth: birth[1][0]) + try: + persistence = _array_handler(persistence) + persistence = _limit_to_max_intervals( + persistence, max_intervals, key=lambda life_time: life_time[1][1] - life_time[1][0] + ) + (min_birth, max_death) = __min_birth_max_death(persistence) + persistence = sorted(persistence, key=lambda birth: birth[1][0]) + except IndexError: + min_birth, max_death = 0.0, 1.0 + pass - (min_birth, max_death) = __min_birth_max_death(persistence) - ind = 0 delta = (max_death - min_birth) * inf_delta # Replace infinity values with max_death + delta for bar code to be more # readable infinity = max_death + delta axis_start = min_birth - delta - # Draw horizontal bars in loop - for interval in reversed(persistence): - if float(interval[1][1]) != float("inf"): - # Finite death case - axes.barh( - ind, - (interval[1][1] - interval[1][0]), - height=0.8, - left=interval[1][0], - alpha=alpha, - color=colormap[interval[0]], - linewidth=0, - ) - else: - # Infinite death case for diagram to be nicer - axes.barh( - ind, - (infinity - interval[1][0]), - height=0.8, - left=interval[1][0], - alpha=alpha, - color=colormap[interval[0]], - linewidth=0, - ) - ind = ind + 1 + + if axes == None: + _, axes = plt.subplots(1, 1) + if colormap == None: + colormap = plt.cm.Set1.colors + + x=[birth for (dim,(birth,death)) in persistence] + 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) if legend: dimensions = list(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", + handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions], loc="lower right", ) axes.set_title("Persistence barcode", fontsize=fontsize) # Ends plot on infinity value and starts a little bit before min_birth - axes.axis([axis_start, infinity, 0, ind]) + if len(x) != 0: + axes.axis([axis_start, infinity, 0, len(x)]) return axes - except ImportError: - print("This function is not available, you may be missing matplotlib.") + except ImportError as import_error: + warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.") def plot_persistence_diagram( @@ -219,14 +218,13 @@ def plot_persistence_diagram( persistence_file="", alpha=0.6, band=0.0, - max_intervals=1000, - max_plots=1000, + max_intervals=1000000, inf_delta=0.1, legend=False, colormap=None, axes=None, fontsize=16, - greyblock=True + greyblock=True, ): """This function plots the persistence diagram from persistence values list, a np.array of shape (N x 2) representing a diagram in a single @@ -244,7 +242,7 @@ def plot_persistence_diagram( :type band: float. :param max_intervals: maximal number of intervals to display. Selected intervals are those with the longest life time. Set it - to 0 to see all. Default value is 1000. + to 0 to see all. Default value is 1000000. :type max_intervals: int. :param inf_delta: Infinity is placed at :code:`((max_death - min_birth) x inf_delta)` above :code:`max_death` value. A reasonable value is @@ -268,47 +266,35 @@ def plot_persistence_diagram( import matplotlib.pyplot as plt import matplotlib.patches as mpatches from matplotlib import rc + if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex(): - plt.rc('text', usetex=True) - plt.rc('font', family='serif') + plt.rc("text", usetex=True) + plt.rc("font", family="serif") else: - plt.rc('text', usetex=False) - plt.rc('font', family='DejaVu Sans') + plt.rc("text", usetex=False) + plt.rc("font", family="DejaVu Sans") if persistence_file != "": if path.isfile(persistence_file): # Reset persistence persistence = [] - diag = read_persistence_intervals_grouped_by_dimension( - persistence_file=persistence_file - ) + diag = read_persistence_intervals_grouped_by_dimension(persistence_file=persistence_file) for key in diag.keys(): for persistence_interval in diag[key]: persistence.append((key, persistence_interval)) else: - print("file " + persistence_file + " not found.") - return None - - persistence = _array_handler(persistence) + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file) - if max_plots != 1000: - print("Deprecated parameter. It has been replaced by max_intervals") - max_intervals = max_plots - - if max_intervals > 0 and max_intervals < len(persistence): - # Sort by life time, then takes only the max_intervals elements - persistence = sorted( - persistence, - key=lambda life_time: life_time[1][1] - life_time[1][0], - reverse=True, - )[:max_intervals] - - if colormap == None: - colormap = plt.cm.Set1.colors - if axes == None: - fig, axes = plt.subplots(1, 1) + try: + persistence = _array_handler(persistence) + persistence = _limit_to_max_intervals( + persistence, max_intervals, key=lambda life_time: life_time[1][1] - life_time[1][0] + ) + min_birth, max_death = __min_birth_max_death(persistence, band) + except IndexError: + min_birth, max_death = 0.0, 1.0 + pass - (min_birth, max_death) = __min_birth_max_death(persistence, band) delta = (max_death - min_birth) * inf_delta # Replace infinity values with max_death + delta for diagram to be more # readable @@ -316,61 +302,56 @@ def plot_persistence_diagram( axis_end = max_death + delta / 2 axis_start = min_birth - delta + if axes == None: + _, axes = plt.subplots(1, 1) + if colormap == None: + colormap = plt.cm.Set1.colors # bootstrap band if band > 0.0: x = np.linspace(axis_start, infinity, 1000) axes.fill_between(x, x, x + band, alpha=alpha, facecolor="red") # lower diag patch if greyblock: - axes.add_patch(mpatches.Polygon([[axis_start, axis_start], [axis_end, axis_start], [axis_end, axis_end]], fill=True, color='lightgrey')) - # Draw points in loop - pts_at_infty = False # Records presence of pts at infty - for interval in reversed(persistence): - if float(interval[1][1]) != float("inf"): - # Finite death case - axes.scatter( - interval[1][0], - interval[1][1], - alpha=alpha, - color=colormap[interval[0]], + axes.add_patch( + mpatches.Polygon( + [[axis_start, axis_start], [axis_end, axis_start], [axis_end, axis_end]], + fill=True, + color="lightgrey", ) - else: - pts_at_infty = True - # Infinite death case for diagram to be nicer - axes.scatter( - interval[1][0], infinity, alpha=alpha, color=colormap[interval[0]] - ) - if pts_at_infty: + ) + # line display of equation : birth = death + axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k") + + x=[birth for (dim,(birth,death)) in persistence] + y=[death if death != float("inf") else infinity for (dim,(birth,death)) in persistence] + c=[colormap[dim] for (dim,(birth,death)) in persistence] + + axes.scatter(x,y,alpha=alpha,color=c) + if float("inf") in (death for (dim,(birth,death)) in persistence): # infinity line and text - axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k") axes.plot([axis_start, axis_end], [infinity, infinity], linewidth=1.0, color="k", alpha=alpha) # Infinity label yt = axes.get_yticks() - yt = yt[np.where(yt < axis_end)] # to avoid ploting ticklabel higher than infinity + yt = yt[np.where(yt < axis_end)] # to avoid ploting ticklabel higher than infinity yt = np.append(yt, infinity) ytl = ["%.3f" % e for e in yt] # to avoid float precision error - ytl[-1] = r'$+\infty$' + ytl[-1] = r"$+\infty$" axes.set_yticks(yt) axes.set_yticklabels(ytl) if legend: dimensions = list(set(item[0] for item in persistence)) - axes.legend( - handles=[ - mpatches.Patch(color=colormap[dim], label=str(dim)) - for dim in dimensions - ] - ) + axes.legend(handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions]) axes.set_xlabel("Birth", fontsize=fontsize) axes.set_ylabel("Death", fontsize=fontsize) axes.set_title("Persistence diagram", fontsize=fontsize) # Ends plot on infinity value and starts a little bit before min_birth - axes.axis([axis_start, axis_end, axis_start, infinity + delta/2]) + axes.axis([axis_start, axis_end, axis_start, infinity + delta / 2]) return axes - except ImportError: - print("This function is not available, you may be missing matplotlib.") + except ImportError as import_error: + warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.") def plot_persistence_density( @@ -384,7 +365,7 @@ def plot_persistence_density( legend=False, axes=None, fontsize=16, - greyblock=False + greyblock=False, ): """This function plots the persistence density from persistence values list, np.array of shape (N x 2) representing a diagram @@ -444,12 +425,13 @@ def plot_persistence_density( import matplotlib.patches as mpatches from scipy.stats import kde from matplotlib import rc + if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex(): - plt.rc('text', usetex=True) - plt.rc('font', family='serif') + plt.rc("text", usetex=True) + plt.rc("font", family="serif") else: - plt.rc('text', usetex=False) - plt.rc('font', family='DejaVu Sans') + plt.rc("text", usetex=False) + plt.rc("font", family="DejaVu Sans") if persistence_file != "": if dimension is None: @@ -460,10 +442,16 @@ def plot_persistence_density( persistence_file=persistence_file, only_this_dim=dimension ) else: - print("file " + persistence_file + " not found.") - return None + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file) + + # default cmap value cannot be done at argument definition level as matplotlib is not yet defined. + if cmap is None: + cmap = plt.cm.hot_r + if axes == None: + _, axes = plt.subplots(1, 1) - if len(persistence) > 0: + try: + # if not read from file but given by an argument persistence = _array_handler(persistence) persistence_dim = np.array( [ @@ -472,47 +460,54 @@ def plot_persistence_density( if (dim_interval[0] == dimension) or (dimension is None) ] ) - - persistence_dim = persistence_dim[np.isfinite(persistence_dim[:, 1])] - if max_intervals > 0 and max_intervals < len(persistence_dim): - # Sort by life time, then takes only the max_intervals elements + persistence_dim = persistence_dim[np.isfinite(persistence_dim[:, 1])] persistence_dim = np.array( - sorted( - persistence_dim, - key=lambda life_time: life_time[1] - life_time[0], - reverse=True, - )[:max_intervals] + _limit_to_max_intervals( + persistence_dim, max_intervals, key=lambda life_time: life_time[1] - life_time[0] + ) ) - # Set as numpy array birth and death (remove undefined values - inf and NaN) - birth = persistence_dim[:, 0] - death = persistence_dim[:, 1] - - # default cmap value cannot be done at argument definition level as matplotlib is not yet defined. - if cmap is None: - cmap = plt.cm.hot_r - if axes == None: - fig, axes = plt.subplots(1, 1) + # Set as numpy array birth and death (remove undefined values - inf and NaN) + birth = persistence_dim[:, 0] + death = persistence_dim[:, 1] + birth_min = birth.min() + birth_max = birth.max() + death_min = death.min() + death_max = death.max() + + # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents + k = kde.gaussian_kde([birth, death], bw_method=bw_method) + xi, yi = np.mgrid[ + birth_min : birth_max : nbins * 1j, death_min : death_max : nbins * 1j, + ] + zi = k(np.vstack([xi.flatten(), yi.flatten()])) + # Make the plot + img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap, shading="auto") + plot_success = True + + # IndexError on empty diagrams, ValueError on only inf death values + except (IndexError, ValueError): + birth_min = 0.0 + birth_max = 1.0 + death_min = 0.0 + death_max = 1.0 + plot_success = False + pass # line display of equation : birth = death - x = np.linspace(death.min(), birth.max(), 1000) + x = np.linspace(death_min, birth_max, 1000) axes.plot(x, x, color="k", linewidth=1.0) - # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents - k = kde.gaussian_kde([birth, death], bw_method=bw_method) - xi, yi = np.mgrid[ - birth.min() : birth.max() : nbins * 1j, - death.min() : death.max() : nbins * 1j, - ] - zi = k(np.vstack([xi.flatten(), yi.flatten()])) - - # Make the plot - img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap) - if greyblock: - axes.add_patch(mpatches.Polygon([[birth.min(), birth.min()], [death.max(), birth.min()], [death.max(), death.max()]], fill=True, color='lightgrey')) + axes.add_patch( + mpatches.Polygon( + [[birth_min, birth_min], [death_max, birth_min], [death_max, death_max]], + fill=True, + color="lightgrey", + ) + ) - if legend: + if plot_success and legend: plt.colorbar(img, ax=axes) axes.set_xlabel("Birth", fontsize=fontsize) @@ -521,7 +516,5 @@ def plot_persistence_density( return axes - except ImportError: - print( - "This function is not available, you may be missing matplotlib and/or scipy." - ) + except ImportError as import_error: + warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.") diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 006a24ed..4f229663 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -45,6 +45,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_interface_full_featured "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>": Simplex_tree_interface_full_featured() nogil + Simplex_tree_interface_full_featured(Simplex_tree_interface_full_featured&) nogil double simplex_filtration(vector[int] simplex) nogil void assign_simplex_filtration(vector[int] simplex, double filtration) nogil void initialize_filtration() nogil @@ -65,6 +66,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]] dgm, double min_persistence) nogil Simplex_tree_interface_full_featured* collapse_edges(int nb_collapse_iteration) nogil except + void reset_filtration(double filtration, int dimension) nogil + bint operator==(Simplex_tree_interface_full_featured) nogil # Iterators over Simplex tree pair[vector[int], double] get_simplex_and_filtration(Simplex_tree_simplex_handle f_simplex) nogil Simplex_tree_simplices_iterator get_simplices_iterator_begin() nogil @@ -74,6 +76,9 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": Simplex_tree_skeleton_iterator get_skeleton_iterator_begin(int dimension) nogil Simplex_tree_skeleton_iterator get_skeleton_iterator_end(int dimension) nogil pair[Simplex_tree_boundary_iterator, Simplex_tree_boundary_iterator] get_boundary_iterators(vector[int] simplex) nogil except + + # Expansion with blockers + ctypedef bool (*blocker_func_t)(vector[int], void *user_data) + void expansion_with_blockers_callback(int dimension, blocker_func_t user_func, void *user_data) cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>": diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index c3720936..2c53a872 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -16,6 +16,9 @@ __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" __license__ = "MIT" +cdef bool callback(vector[int] simplex, void *blocker_func): + return (<object>blocker_func)(simplex) + # SimplexTree python interface cdef class SimplexTree: """The simplex tree is an efficient and flexible data structure for @@ -38,13 +41,29 @@ cdef class SimplexTree: cdef Simplex_tree_persistence_interface * pcohptr # Fake constructor that does nothing but documenting the constructor - def __init__(self): + def __init__(self, other = None): """SimplexTree constructor. + + :param other: If `other` is `None` (default value), an empty `SimplexTree` is created. + If `other` is a `SimplexTree`, the `SimplexTree` is constructed from a deep copy of `other`. + :type other: SimplexTree (Optional) + :returns: An empty or a copy simplex tree. + :rtype: SimplexTree + + :raises TypeError: In case `other` is neither `None`, nor a `SimplexTree`. + :note: If the `SimplexTree` is a copy, the persistence information is not copied. If you need it in the clone, + you have to call :func:`compute_persistence` on it even if you had already computed it in the original. """ # The real cython constructor - def __cinit__(self): - self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured()) + def __cinit__(self, other = None): + if other: + if isinstance(other, SimplexTree): + self.thisptr = _get_copy_intptr(other) + else: + raise TypeError("`other` argument requires to be of type `SimplexTree`, or `None`.") + else: + self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured()) def __dealloc__(self): cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr() @@ -63,6 +82,21 @@ cdef class SimplexTree: """ return self.pcohptr != NULL + def copy(self): + """ + :returns: A simplex tree that is a deep copy of itself. + :rtype: SimplexTree + + :note: The persistence information is not copied. If you need it in the clone, you have to call + :func:`compute_persistence` on it even if you had already computed it in the original. + """ + stree = SimplexTree() + stree.thisptr = _get_copy_intptr(self) + return stree + + def __deepcopy__(self): + return self.copy() + def filtration(self, simplex): """This function returns the filtration value for a given N-simplex in this simplicial complex, or +infinity if it is not in the complex. @@ -443,6 +477,27 @@ cdef class SimplexTree: persistence_result = self.pcohptr.get_persistence() return self.get_ptr().compute_extended_persistence_subdiagrams(persistence_result, min_persistence) + def expansion_with_blocker(self, max_dim, blocker_func): + """Expands the Simplex_tree containing only a graph. Simplices corresponding to cliques in the graph are added + incrementally, faces before cofaces, unless the simplex has dimension larger than `max_dim` or `blocker_func` + returns `True` for this simplex. + + The function identifies a candidate simplex whose faces are all already in the complex, inserts it with a + filtration value corresponding to the maximum of the filtration values of the faces, then calls `blocker_func` + with this new simplex (represented as a list of int). If `blocker_func` returns `True`, the simplex is removed, + otherwise it is kept. The algorithm then proceeds with the next candidate. + + .. warning:: + Several candidates of the same dimension may be inserted simultaneously before calling `block_simplex`, so + if you examine the complex in `block_simplex`, you may hit a few simplices of the same dimension that have + not been vetted by `block_simplex` yet, or have already been rejected but not yet removed. + + :param max_dim: Expansion maximal dimension value. + :type max_dim: int + :param blocker_func: Blocker oracle. + :type blocker_func: Callable[[List[int]], bool] + """ + self.get_ptr().expansion_with_blockers_callback(max_dim, callback, <void*>blocker_func) def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False): """This function computes and returns the persistence of the simplicial complex. @@ -621,18 +676,17 @@ cdef class SimplexTree: return (normal0, normals, infinite0, infinites) def collapse_edges(self, nb_iterations = 1): - """Assuming the simplex tree is a 1-skeleton graph, this method collapse edges (simplices of higher dimension - are ignored) and resets the simplex tree from the remaining edges. - A good candidate is to build a simplex tree on top of a :class:`~gudhi.RipsComplex` of dimension 1 before - collapsing edges + """Assuming the complex is a graph (simplices of higher dimension are ignored), this method implicitly + interprets it as the 1-skeleton of a flag complex, and replaces it with another (smaller) graph whose + expansion has the same persistent homology, using a technique known as edge collapses + (see :cite:`edgecollapsearxiv`). + + A natural application is to get a simplex tree of dimension 1 from :class:`~gudhi.RipsComplex`, + then collapse edges, perform :meth:`expansion()` and finally compute persistence (cf. :download:`rips_complex_edge_collapse_example.py <../example/rips_complex_edge_collapse_example.py>`). - For implementation details, please refer to :cite:`edgecollapsesocg2020`. :param nb_iterations: The number of edge collapse iterations to perform. Default is 1. :type nb_iterations: int - - :note: collapse_edges method requires `Eigen <installation.html#eigen>`_ >= 3.1.0 and an exception is thrown - if this method is not available. """ # Backup old pointer cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr() @@ -642,3 +696,13 @@ cdef class SimplexTree: self.thisptr = <intptr_t>(ptr.collapse_edges(nb_iter)) # Delete old pointer del ptr + + def __eq__(self, other:SimplexTree): + """Test for structural equality + :returns: True if the 2 simplex trees are equal, False otherwise. + :rtype: bool + """ + return dereference(self.get_ptr()) == dereference(other.get_ptr()) + +cdef intptr_t _get_copy_intptr(SimplexTree stree) nogil: + return <intptr_t>(new Simplex_tree_interface_full_featured(dereference(stree.get_ptr()))) diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py index 0541572b..16f63c3d 100644 --- a/src/python/gudhi/weighted_rips_complex.py +++ b/src/python/gudhi/weighted_rips_complex.py @@ -12,9 +12,11 @@ from gudhi import SimplexTree class WeightedRipsComplex: """ Class to generate a weighted Rips complex from a distance matrix and weights on vertices, - in the way described in :cite:`dtmfiltrations`. + in the way described in :cite:`dtmfiltrations` with `p=1`. The filtration value of vertex `i` is `2*weights[i]`, + and the filtration value of edge `ij` is `distance_matrix[i][j]+weights[i]+weights[j]`, + or the maximum of the filtrations of its extremities, whichever is largest. Remark that all the filtration values are doubled compared to the definition in the paper - for the consistency with RipsComplex. + for consistency with RipsComplex. """ def __init__(self, distance_matrix, diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 671af4a4..469b91ce 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -57,6 +57,16 @@ class Alpha_complex_interface { alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); } + static void set_float_relative_precision(double precision) { + // cf. Exact_alpha_complex_dD kernel type in Alpha_complex_factory.h + CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>::FT::set_relative_precision_of_to_double(precision); + } + + static double get_float_relative_precision() { + // cf. Exact_alpha_complex_dD kernel type in Alpha_complex_factory.h + return CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>::FT::get_relative_precision_of_to_double(); + } + private: std::unique_ptr<Abstract_alpha_complex> alpha_ptr_; }; diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 629f6083..7f9b0067 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -15,9 +15,7 @@ #include <gudhi/distance_functions.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Points_off_io.h> -#ifdef GUDHI_USE_EIGEN3 #include <gudhi/Flag_complex_edge_collapser.h> -#endif #include <iostream> #include <vector> @@ -42,6 +40,7 @@ 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; + typedef bool (*blocker_func_t)(Simplex simplex, void *user_data); public: @@ -164,7 +163,6 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { } Simplex_tree_interface* collapse_edges(int nb_collapse_iteration) { -#ifdef GUDHI_USE_EIGEN3 using Filtered_edge = std::tuple<Vertex_handle, Vertex_handle, Filtration_value>; std::vector<Filtered_edge> edges; for (Simplex_handle sh : Base::skeleton_simplex_range(1)) { @@ -178,7 +176,7 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { } for (int iteration = 0; iteration < nb_collapse_iteration; iteration++) { - edges = Gudhi::collapse::flag_complex_collapse_edges(edges); + edges = Gudhi::collapse::flag_complex_collapse_edges(std::move(edges)); } Simplex_tree_interface* collapsed_stree_ptr = new Simplex_tree_interface(); // Copy the original 0-skeleton @@ -190,9 +188,13 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { collapsed_stree_ptr->insert({std::get<0>(remaining_edge), std::get<1>(remaining_edge)}, std::get<2>(remaining_edge)); } return collapsed_stree_ptr; -#else - throw std::runtime_error("Unable to collapse edges as it requires Eigen3 >= 3.1.0."); -#endif + } + + void expansion_with_blockers_callback(int dimension, blocker_func_t user_func, void *user_data) { + Base::expansion_with_blockers(dimension, [&](Simplex_handle sh){ + Simplex simplex(Base::simplex_vertex_range(sh).begin(), Base::simplex_vertex_range(sh).end()); + return user_func(simplex, user_data); + }); } // Iterator over the simplex tree diff --git a/src/python/pyproject.toml b/src/python/pyproject.toml index a9fb4985..55b64466 100644 --- a/src/python/pyproject.toml +++ b/src/python/pyproject.toml @@ -1,3 +1,3 @@ [build-system] -requires = ["setuptools", "wheel", "numpy>=1.15.0", "cython", "pybind11"] +requires = ["setuptools>=24.2.0", "wheel", "numpy>=1.15.0", "cython>=0.27", "pybind11"] build-backend = "setuptools.build_meta" diff --git a/src/python/setup.py.in b/src/python/setup.py.in index 23746998..2c67c2c5 100644 --- a/src/python/setup.py.in +++ b/src/python/setup.py.in @@ -5,6 +5,7 @@ Copyright (C) 2019 Inria Modification(s): + - 2021/12 Vincent Rouvreau: Python 3.5 as minimal version - YYYY/MM Author: Description of the modification """ @@ -43,7 +44,7 @@ for module in cython_modules: include_dirs=include_dirs, runtime_library_dirs=runtime_library_dirs,)) -ext_modules = cythonize(ext_modules, compiler_directives={'language_level': str(sys.version_info[0])}) +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)] @@ -86,6 +87,7 @@ setup( long_description_content_type='text/x-rst', long_description=long_description, ext_modules = ext_modules, + python_requires='>=3.5.0', install_requires = ['numpy >= 1.15.0',], package_data={"": ["*.dll"], }, ) diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py index f15284f3..f81e6137 100755 --- a/src/python/test/test_alpha_complex.py +++ b/src/python/test/test_alpha_complex.py @@ -286,3 +286,30 @@ def _weighted_doc_example(precision): def test_weighted_doc_example(): for precision in ['fast', 'safe', 'exact']: _weighted_doc_example(precision) + +def test_float_relative_precision(): + assert AlphaComplex.get_float_relative_precision() == 1e-5 + # Must be > 0. + with pytest.raises(ValueError): + AlphaComplex.set_float_relative_precision(0.) + # Must be < 1. + with pytest.raises(ValueError): + AlphaComplex.set_float_relative_precision(1.) + + points = [[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]] + st = AlphaComplex(points=points).create_simplex_tree() + filtrations = list(st.get_filtration()) + + # Get a better precision + AlphaComplex.set_float_relative_precision(1e-15) + assert AlphaComplex.get_float_relative_precision() == 1e-15 + + st = AlphaComplex(points=points).create_simplex_tree() + filtrations_better_resolution = list(st.get_filtration()) + + assert len(filtrations) == len(filtrations_better_resolution) + for idx in range(len(filtrations)): + # check simplex is the same + assert filtrations[idx][0] == filtrations_better_resolution[idx][0] + # check filtration is about the same with a relative precision of the worst case + assert filtrations[idx][1] == pytest.approx(filtrations_better_resolution[idx][1], rel=1e-5) diff --git a/src/python/test/test_persistence_graphical_tools.py b/src/python/test/test_persistence_graphical_tools.py new file mode 100644 index 00000000..c19836b7 --- /dev/null +++ b/src/python/test/test_persistence_graphical_tools.py @@ -0,0 +1,121 @@ +""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + Author(s): Vincent Rouvreau + + Copyright (C) 2021 Inria + + Modification(s): + - YYYY/MM Author: Description of the modification +""" + +import gudhi as gd +import numpy as np +import matplotlib as plt +import pytest + + +def test_array_handler(): + diags = np.array([[1, 2], [3, 4], [5, 6]], float) + arr_diags = gd.persistence_graphical_tools._array_handler(diags) + for idx in range(len(diags)): + assert arr_diags[idx][0] == 0 + np.testing.assert_array_equal(arr_diags[idx][1], diags[idx]) + + diags = [(1.0, 2.0), (3.0, 4.0), (5.0, 6.0)] + arr_diags = gd.persistence_graphical_tools._array_handler(diags) + for idx in range(len(diags)): + assert arr_diags[idx][0] == 0 + assert arr_diags[idx][1] == diags[idx] + + diags = [(0, (1.0, 2.0)), (0, (3.0, 4.0)), (0, (5.0, 6.0))] + assert gd.persistence_graphical_tools._array_handler(diags) == diags + + +def test_min_birth_max_death(): + diags = [ + (0, (0.0, float("inf"))), + (0, (0.0983494, float("inf"))), + (0, (0.0, 0.122545)), + (0, (0.0, 0.12047)), + (0, (0.0, 0.118398)), + (0, (0.118398, 1.0)), + (0, (0.0, 0.117908)), + (0, (0.0, 0.112307)), + (0, (0.0, 0.107535)), + (0, (0.0, 0.106382)), + ] + assert gd.persistence_graphical_tools.__min_birth_max_death(diags) == (0.0, 1.0) + assert gd.persistence_graphical_tools.__min_birth_max_death(diags, band=4.0) == (0.0, 5.0) + + +def test_limit_min_birth_max_death(): + diags = [ + (0, (2.0, float("inf"))), + (0, (2.0, float("inf"))), + ] + assert gd.persistence_graphical_tools.__min_birth_max_death(diags) == (2.0, 3.0) + assert gd.persistence_graphical_tools.__min_birth_max_death(diags, band=4.0) == (2.0, 6.0) + + +def test_limit_to_max_intervals(): + diags = [ + (0, (0.0, float("inf"))), + (0, (0.0983494, float("inf"))), + (0, (0.0, 0.122545)), + (0, (0.0, 0.12047)), + (0, (0.0, 0.118398)), + (0, (0.118398, 1.0)), + (0, (0.0, 0.117908)), + (0, (0.0, 0.112307)), + (0, (0.0, 0.107535)), + (0, (0.0, 0.106382)), + ] + # check no warnings if max_intervals equals to the diagrams number + with pytest.warns(None) as record: + 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: + truncated_diags = gd.persistence_graphical_tools._limit_to_max_intervals( + diags, 5, key=lambda life_time: life_time[1][1] - life_time[1][0] + ) + # check diagrams are truncated and sorted by life time + assert truncated_diags == [ + (0, (0.0, float("inf"))), + (0, (0.0983494, float("inf"))), + (0, (0.118398, 1.0)), + (0, (0.0, 0.122545)), + (0, (0.0, 0.12047)), + ] + assert len(record) == 1 + + +def _limit_plot_persistence(function): + pplot = function(persistence=[]) + assert isinstance(pplot, plt.axes.SubplotBase) + pplot = function(persistence=[], legend=True) + assert isinstance(pplot, plt.axes.SubplotBase) + pplot = function(persistence=[(0, float("inf"))]) + assert isinstance(pplot, plt.axes.SubplotBase) + pplot = function(persistence=[(0, float("inf"))], legend=True) + assert isinstance(pplot, plt.axes.SubplotBase) + + +def test_limit_plot_persistence(): + for function in [gd.plot_persistence_barcode, gd.plot_persistence_diagram, gd.plot_persistence_density]: + _limit_plot_persistence(function) + + +def _non_existing_persistence_file(function): + with pytest.raises(FileNotFoundError): + function(persistence_file="pouetpouettralala.toubiloubabdou") + + +def test_non_existing_persistence_file(): + for function in [gd.plot_persistence_barcode, gd.plot_persistence_diagram, gd.plot_persistence_density]: + _non_existing_persistence_file(function) diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 31c46213..688f4fd6 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -8,7 +8,7 @@ - YYYY/MM Author: Description of the modification """ -from gudhi import SimplexTree, __GUDHI_USE_EIGEN3 +from gudhi import SimplexTree import numpy as np import pytest @@ -354,16 +354,11 @@ def test_collapse_edges(): assert st.num_simplices() == 10 - if __GUDHI_USE_EIGEN3: - st.collapse_edges() - assert st.num_simplices() == 9 - assert st.find([1, 3]) == False - for simplex in st.get_skeleton(0): - assert simplex[1] == 1. - else: - # If no Eigen3, collapse_edges throws an exception - with pytest.raises(RuntimeError): - st.collapse_edges() + st.collapse_edges() + assert st.num_simplices() == 9 + assert st.find([0, 2]) == False # [1, 3] would be fine as well + for simplex in st.get_skeleton(0): + assert simplex[1] == 1. def test_reset_filtration(): st = SimplexTree() @@ -447,4 +442,116 @@ def test_persistence_intervals_in_dimension(): assert np.array_equal(H2, np.array([[ 0., float("inf")]])) # Test empty case assert st.persistence_intervals_in_dimension(3).shape == (0, 2) -
\ No newline at end of file + +def test_equality_operator(): + st1 = SimplexTree() + st2 = SimplexTree() + + assert st1 == st2 + + st1.insert([1,2,3], 4.) + assert st1 != st2 + + st2.insert([1,2,3], 4.) + assert st1 == st2 + +def test_simplex_tree_deep_copy(): + st = SimplexTree() + st.insert([1, 2, 3], 0.) + # compute persistence only on the original + st.compute_persistence() + + st_copy = st.copy() + assert st_copy == st + st_filt_list = list(st.get_filtration()) + + # check persistence is not copied + assert st.__is_persistence_defined() == True + assert st_copy.__is_persistence_defined() == False + + # remove something in the copy and check the copy is included in the original + st_copy.remove_maximal_simplex([1, 2, 3]) + a_filt_list = list(st_copy.get_filtration()) + assert len(a_filt_list) < len(st_filt_list) + + for a_splx in a_filt_list: + assert a_splx in st_filt_list + + # test double free + del st + del st_copy + +def test_simplex_tree_deep_copy_constructor(): + st = SimplexTree() + st.insert([1, 2, 3], 0.) + # compute persistence only on the original + st.compute_persistence() + + st_copy = SimplexTree(st) + assert st_copy == st + st_filt_list = list(st.get_filtration()) + + # check persistence is not copied + assert st.__is_persistence_defined() == True + assert st_copy.__is_persistence_defined() == False + + # remove something in the copy and check the copy is included in the original + st_copy.remove_maximal_simplex([1, 2, 3]) + a_filt_list = list(st_copy.get_filtration()) + assert len(a_filt_list) < len(st_filt_list) + + for a_splx in a_filt_list: + assert a_splx in st_filt_list + + # test double free + del st + del st_copy + +def test_simplex_tree_constructor_exception(): + with pytest.raises(TypeError): + st = SimplexTree(other = "Construction from a string shall raise an exception") + +def test_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) + + def blocker(simplex): + try: + # Block all simplices that countains vertex 6 + simplex.index(6) + print(simplex, ' is blocked') + return True + except ValueError: + print(simplex, ' is accepted') + st.assign_filtration(simplex, st.filtration(simplex) + 1.) + 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. + + 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. |