diff options
Diffstat (limited to 'src/python')
| -rw-r--r-- | src/python/CMakeLists.txt | 14 | ||||
| -rw-r--r-- | src/python/doc/alpha_complex_ref.rst | 7 | ||||
| -rw-r--r-- | src/python/doc/alpha_complex_sum.inc | 24 | ||||
| -rw-r--r-- | src/python/doc/alpha_complex_user.rst | 110 | ||||
| -rwxr-xr-x | src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py | 55 | ||||
| -rwxr-xr-x | src/python/example/alpha_rips_persistence_bottleneck_distance.py | 110 | ||||
| -rwxr-xr-x | src/python/example/plot_alpha_complex.py | 5 | ||||
| -rw-r--r-- | src/python/gudhi/alpha_complex.pyx | 49 | ||||
| -rw-r--r-- | src/python/gudhi/alpha_complex_3d.pyx | 129 | ||||
| -rw-r--r-- | src/python/include/Alpha_complex_factory.h | 111 | ||||
| -rw-r--r-- | src/python/include/Alpha_complex_interface.h | 54 | ||||
| -rw-r--r-- | src/python/include/Alpha_complex_interface_3d.h | 71 | ||||
| -rwxr-xr-x | src/python/test/test_alpha_complex.py | 87 | ||||
| -rwxr-xr-x | src/python/test/test_reader_utils.py | 33 |
14 files changed, 615 insertions, 244 deletions
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 73303a24..0739d7b5 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -61,6 +61,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'subsampling', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'tangential_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'alpha_complex', ") + set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'alpha_complex_3d', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'euclidean_witness_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'euclidean_strong_witness_complex', ") # Modules that should not be auto-imported in __init__.py @@ -155,12 +156,15 @@ if(PYTHONINTERP_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'nerve_gic', ") endif () if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) - set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") + set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex_3d', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'subsampling', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'tangential_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_witness_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_strong_witness_complex', ") endif () + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") + endif () if(CGAL_FOUND) # Add CGAL compilation args @@ -343,13 +347,15 @@ if(PYTHONINTERP_FOUND) # Test examples - if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) # Bottleneck and Alpha add_test(NAME alpha_rips_persistence_bottleneck_distance_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/alpha_rips_persistence_bottleneck_distance.py" -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -t 0.15 -d 3) + endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) # Tangential add_test(NAME tangential_complex_plain_homology_from_off_file_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} @@ -419,7 +425,7 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_cover_complex) endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) - if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) # Alpha add_test(NAME alpha_complex_from_points_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} @@ -431,7 +437,7 @@ if(PYTHONINTERP_FOUND) ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/alpha_complex_diagram_persistence_from_off_file_example.py" --no-diagram -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off) add_gudhi_py_test(test_alpha_complex) - endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) # Euclidean witness diff --git a/src/python/doc/alpha_complex_ref.rst b/src/python/doc/alpha_complex_ref.rst index 7da79543..49321368 100644 --- a/src/python/doc/alpha_complex_ref.rst +++ b/src/python/doc/alpha_complex_ref.rst @@ -9,6 +9,11 @@ Alpha complex reference manual .. autoclass:: gudhi.AlphaComplex :members: :undoc-members: - :show-inheritance: .. automethod:: gudhi.AlphaComplex.__init__ + +.. autoclass:: gudhi.AlphaComplex3D + :members: + :undoc-members: + + .. automethod:: gudhi.AlphaComplex3D.__init__ diff --git a/src/python/doc/alpha_complex_sum.inc b/src/python/doc/alpha_complex_sum.inc index aeab493f..5c76fd54 100644 --- a/src/python/doc/alpha_complex_sum.inc +++ b/src/python/doc/alpha_complex_sum.inc @@ -1,15 +1,15 @@ .. table:: :widths: 30 40 30 - +----------------------------------------------------------------+-------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------------+ - | .. figure:: | Alpha complex is a simplicial complex constructed from the finite | :Author: Vincent Rouvreau | - | ../../doc/Alpha_complex/alpha_complex_representation.png | cells of a Delaunay Triangulation. It has the same persistent homology | | - | :alt: Alpha complex representation | as the Čech complex and is significantly smaller. | :Since: GUDHI 2.0.0 | - | :figclass: align-center | | | - | | | :License: MIT (`GPL v3 </licensing/>`_) | - | | | | - | | | :Requires: `Eigen <installation.html#eigen>`_ :math:`\geq` 3.1.0 and `CGAL <installation.html#cgal>`_ :math:`\geq` 4.11.0 | - | | | | - +----------------------------------------------------------------+-------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------------+ - | * :doc:`alpha_complex_user` | * :doc:`alpha_complex_ref` | - +----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ + +----------------------------------------------------------------+-------------------------------------------------------------------------+------------------------------------------------------------------------------------------------------------------------+ + | .. figure:: | Alpha complex is a simplicial complex constructed from the finite | :Author: Vincent Rouvreau | + | ../../doc/Alpha_complex/alpha_complex_representation.png | cells of a Delaunay Triangulation. It has the same persistent homology | | + | :alt: Alpha complex representation | as the Čech complex and is significantly smaller. | :Since: GUDHI 2.0.0 | + | :figclass: align-center | | | + | | | :License: MIT (`GPL v3 </licensing/>`_) | + | | | | + | | | :Requires: `Eigen <installation.html#eigen>`_ :math:`\geq` 3.1.0 and `CGAL <installation.html#cgal>`_ :math:`\geq` 5.1 | + | | | | + +----------------------------------------------------------------+-------------------------------------------------------------------------+------------------------------------------------------------------------------------------------------------------------+ + | * :doc:`alpha_complex_user` | * :doc:`alpha_complex_ref` | + +----------------------------------------------------------------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst index fffcb3db..13187f1f 100644 --- a/src/python/doc/alpha_complex_user.rst +++ b/src/python/doc/alpha_complex_user.rst @@ -9,7 +9,7 @@ Definition .. include:: alpha_complex_sum.inc -:doc:`AlphaComplex <alpha_complex_ref>` is constructing a :doc:`SimplexTree <simplex_tree_ref>` using +:class:`~gudhi.AlphaComplex` is constructing a :doc:`SimplexTree <simplex_tree_ref>` using `Delaunay Triangulation <http://doc.cgal.org/latest/Triangulation/index.html#Chapter_Triangulations>`_ :cite:`cgal:hdj-t-19b` from the `Computational Geometry Algorithms Library <http://www.cgal.org/>`_ :cite:`cgal:eb-19b`. @@ -44,23 +44,22 @@ This example builds the alpha-complex from the given points: .. testcode:: - import gudhi - alpha_complex = gudhi.AlphaComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]]) + from gudhi import AlphaComplex + ac = AlphaComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]]) + + stree = ac.create_simplex_tree() + print('Alpha complex is of dimension ', stree.dimension(), ' - ', + stree.num_simplices(), ' simplices - ', stree.num_vertices(), ' vertices.') - simplex_tree = alpha_complex.create_simplex_tree() - result_str = 'Alpha complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ - repr(simplex_tree.num_simplices()) + ' simplices - ' + \ - repr(simplex_tree.num_vertices()) + ' vertices.' - print(result_str) fmt = '%s -> %.2f' - for filtered_value in simplex_tree.get_filtration(): + for filtered_value in stree.get_filtration(): print(fmt % tuple(filtered_value)) The output is: .. testoutput:: - Alpha complex is of dimension 2 - 25 simplices - 7 vertices. + Alpha complex is of dimension 2 - 25 simplices - 7 vertices. [0] -> 0.00 [1] -> 0.00 [2] -> 0.00 @@ -174,11 +173,76 @@ of speed-up, since we always first build the full filtered complex, so it is rec :paramref:`~gudhi.AlphaComplex.create_simplex_tree.max_alpha_square`. In the following example, a threshold of :math:`\alpha^2 = 32.0` is used. +Weighted specific 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>`_ +of Delaunay. + +This example builds the CGAL weighted alpha shapes from a small molecule, and initializes the alpha complex with +it. This example is taken from +`CGAL 3d weighted alpha shapes <https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title13>`_. + +Then, it is asked to display information about the alpha complex. + +.. testcode:: + + from gudhi import AlphaComplex + wgt_ac = AlphaComplex(points=[[ 1., -1., -1.], + [-1., 1., -1.], + [-1., -1., 1.], + [ 1., 1., 1.], + [ 2., 2., 2.]], + weights = [4., 4., 4., 4., 1.]) + + stree = wgt_ac.create_simplex_tree() + print('Weighted alpha complex is of dimension ', stree.dimension(), ' - ', + stree.num_simplices(), ' simplices - ', stree.num_vertices(), ' vertices.') + fmt = '%s -> %.2f' + for filtered_value in stree.get_filtration(): + print(fmt % tuple(filtered_value)) + +The output is: + +.. testoutput:: + + Weighted alpha complex is of dimension 3 - 29 simplices - 5 vertices. + [0] -> -4.00 + [1] -> -4.00 + [2] -> -4.00 + [3] -> -4.00 + [0, 1] -> -2.00 + [0, 2] -> -2.00 + [1, 2] -> -2.00 + [0, 3] -> -2.00 + [1, 3] -> -2.00 + [2, 3] -> -2.00 + [0, 2, 3] -> -1.33 + [1, 2, 3] -> -1.33 + [0, 1, 2] -> -1.33 + [0, 1, 3] -> -1.33 + [0, 1, 2, 3] -> -1.00 + [4] -> -1.00 + [3, 4] -> -1.00 + [0, 4] -> 23.00 + [1, 4] -> 23.00 + [2, 4] -> 23.00 + [0, 3, 4] -> 23.00 + [1, 3, 4] -> 23.00 + [2, 3, 4] -> 23.00 + [0, 1, 4] -> 95.00 + [0, 2, 4] -> 95.00 + [1, 2, 4] -> 95.00 + [0, 1, 3, 4] -> 95.00 + [0, 2, 3, 4] -> 95.00 + [1, 2, 3, 4] -> 95.00 Example from OFF file ^^^^^^^^^^^^^^^^^^^^^ -This example builds the alpha complex from 300 random points on a 2-torus. +This example builds the alpha complex from 300 random points on a 2-torus, given by an +`OFF file <fileformats.html#off-file-format>`_. Then, it computes the persistence diagram and displays it: @@ -186,14 +250,18 @@ Then, it computes the persistence diagram and displays it: :include-source: import matplotlib.pyplot as plt - import gudhi - alpha_complex = gudhi.AlphaComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/tore3D_300.off') - simplex_tree = alpha_complex.create_simplex_tree() - result_str = 'Alpha complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ - repr(simplex_tree.num_simplices()) + ' simplices - ' + \ - repr(simplex_tree.num_vertices()) + ' vertices.' - print(result_str) - diag = simplex_tree.persistence() - gudhi.plot_persistence_diagram(diag) + import gudhi as gd + off_file = gd.__root_source_dir__ + '/data/points/tore3D_300.off' + points = gd.read_points_from_off_file(off_file = off_file) + stree = gd.AlphaComplex(points = points).create_simplex_tree() + dgm = stree.persistence() + gd.plot_persistence_diagram(dgm, legend = True) plt.show() + +3d specific version +^^^^^^^^^^^^^^^^^^^ + +:Requires: `Eigen <installation.html#eigen>`_ :math:`\geq` 3.1.0 and `CGAL <installation.html#cgal>`_ :math:`\geq` 4.11.0. + +A specific module for Alpha complex is available in 3d (cf. :class:`~gudhi.AlphaComplex3D`) and +allows to construct standard and weighted versions of alpha complexes.
\ No newline at end of file diff --git a/src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py b/src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py index fe03be31..c96121a6 100755 --- a/src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py +++ b/src/python/example/alpha_complex_diagram_persistence_from_off_file_example.py @@ -1,9 +1,7 @@ #!/usr/bin/env python import argparse -import errno -import os -import gudhi +import gudhi as gd """ This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. @@ -41,33 +39,24 @@ parser.add_argument( args = parser.parse_args() -with open(args.file, "r") as f: - first_line = f.readline() - if (first_line == "OFF\n") or (first_line == "nOFF\n"): - print("##############################################################") - print("AlphaComplex creation from points read in a OFF file") - - alpha_complex = gudhi.AlphaComplex(off_file=args.file) - if args.max_alpha_square is not None: - print("with max_edge_length=", args.max_alpha_square) - simplex_tree = alpha_complex.create_simplex_tree( - max_alpha_square=args.max_alpha_square - ) - else: - simplex_tree = alpha_complex.create_simplex_tree() - - print("Number of simplices=", simplex_tree.num_simplices()) - - diag = simplex_tree.persistence() - - print("betti_numbers()=", simplex_tree.betti_numbers()) - - if args.no_diagram == False: - import matplotlib.pyplot as plot - gudhi.plot_persistence_diagram(diag, band=args.band) - plot.show() - else: - raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), - args.file) - - f.close() +print("##############################################################") +print("AlphaComplex creation from points read in a OFF file") + +points = gd.read_points_from_off_file(off_file = args.file) +alpha_complex = gd.AlphaComplex(points = points) +if args.max_alpha_square is not None: + print("with max_edge_length=", args.max_alpha_square) + simplex_tree = alpha_complex.create_simplex_tree( + max_alpha_square=args.max_alpha_square + ) +else: + simplex_tree = alpha_complex.create_simplex_tree() + +print("Number of simplices=", simplex_tree.num_simplices()) + +diag = simplex_tree.persistence() +print("betti_numbers()=", simplex_tree.betti_numbers()) +if args.no_diagram == False: + import matplotlib.pyplot as plot + gd.plot_persistence_diagram(diag, band=args.band) + plot.show() diff --git a/src/python/example/alpha_rips_persistence_bottleneck_distance.py b/src/python/example/alpha_rips_persistence_bottleneck_distance.py index 3e12b0d5..6b97fb3b 100755 --- a/src/python/example/alpha_rips_persistence_bottleneck_distance.py +++ b/src/python/example/alpha_rips_persistence_bottleneck_distance.py @@ -1,10 +1,8 @@ #!/usr/bin/env python -import gudhi +import gudhi as gd import argparse import math -import errno -import os import numpy as np """ This file is part of the Gudhi Library - https://gudhi.inria.fr/ - @@ -37,70 +35,60 @@ parser.add_argument("-t", "--threshold", type=float, default=0.5) parser.add_argument("-d", "--max_dimension", type=int, default=1) args = parser.parse_args() -with open(args.file, "r") as f: - first_line = f.readline() - if (first_line == "OFF\n") or (first_line == "nOFF\n"): - point_cloud = gudhi.read_points_from_off_file(off_file=args.file) - print("##############################################################") - print("RipsComplex creation from points read in a OFF file") +point_cloud = gd.read_points_from_off_file(off_file=args.file) +print("##############################################################") +print("RipsComplex creation from points read in a OFF file") - message = "RipsComplex with max_edge_length=" + repr(args.threshold) - print(message) +message = "RipsComplex with max_edge_length=" + repr(args.threshold) +print(message) - rips_complex = gudhi.RipsComplex( - points=point_cloud, max_edge_length=args.threshold - ) - - rips_stree = rips_complex.create_simplex_tree( - max_dimension=args.max_dimension) - - message = "Number of simplices=" + repr(rips_stree.num_simplices()) - print(message) - - rips_stree.compute_persistence() - - print("##############################################################") - print("AlphaComplex creation from points read in a OFF file") - - message = "AlphaComplex with max_edge_length=" + repr(args.threshold) - print(message) - - alpha_complex = gudhi.AlphaComplex(points=point_cloud) - alpha_stree = alpha_complex.create_simplex_tree( - max_alpha_square=(args.threshold * args.threshold) - ) - - message = "Number of simplices=" + repr(alpha_stree.num_simplices()) - print(message) +rips_complex = gd.RipsComplex( + points=point_cloud, max_edge_length=args.threshold +) - alpha_stree.compute_persistence() +rips_stree = rips_complex.create_simplex_tree( + max_dimension=args.max_dimension) - max_b_distance = 0.0 - for dim in range(args.max_dimension): - # Alpha persistence values needs to be transform because filtration - # values are alpha square values - alpha_intervals = np.sqrt(alpha_stree.persistence_intervals_in_dimension(dim)) +message = "Number of simplices=" + repr(rips_stree.num_simplices()) +print(message) - rips_intervals = rips_stree.persistence_intervals_in_dimension(dim) - bottleneck_distance = gudhi.bottleneck_distance( - rips_intervals, alpha_intervals - ) - message = ( - "In dimension " - + repr(dim) - + ", bottleneck distance = " - + repr(bottleneck_distance) - ) - print(message) - max_b_distance = max(bottleneck_distance, max_b_distance) +rips_stree.compute_persistence() - print("==============================================================") - message = "Bottleneck distance is " + repr(max_b_distance) - print(message) +print("##############################################################") +print("AlphaComplex creation from points read in a OFF file") - else: - raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), - args.file) +message = "AlphaComplex with max_edge_length=" + repr(args.threshold) +print(message) +alpha_complex = gd.AlphaComplex(points=point_cloud) +alpha_stree = alpha_complex.create_simplex_tree( + max_alpha_square=(args.threshold * args.threshold) +) - f.close() +message = "Number of simplices=" + repr(alpha_stree.num_simplices()) +print(message) + +alpha_stree.compute_persistence() + +max_b_distance = 0.0 +for dim in range(args.max_dimension): + # Alpha persistence values needs to be transform because filtration + # values are alpha square values + alpha_intervals = np.sqrt(alpha_stree.persistence_intervals_in_dimension(dim)) + + rips_intervals = rips_stree.persistence_intervals_in_dimension(dim) + bottleneck_distance = gd.bottleneck_distance( + rips_intervals, alpha_intervals + ) + message = ( + "In dimension " + + repr(dim) + + ", bottleneck distance = " + + repr(bottleneck_distance) + ) + print(message) + max_b_distance = max(bottleneck_distance, max_b_distance) + +print("==============================================================") +message = "Bottleneck distance is " + repr(max_b_distance) +print(message) diff --git a/src/python/example/plot_alpha_complex.py b/src/python/example/plot_alpha_complex.py index 99c18a7c..0924619b 100755 --- a/src/python/example/plot_alpha_complex.py +++ b/src/python/example/plot_alpha_complex.py @@ -1,8 +1,9 @@ #!/usr/bin/env python import numpy as np -import gudhi -ac = gudhi.AlphaComplex(off_file='../../data/points/tore3D_1307.off') +import gudhi as gd +points = gd.read_points_from_off_file(off_file = '../../data/points/tore3D_1307.off') +ac = gd.AlphaComplex(points = points) st = ac.create_simplex_tree() points = np.array([ac.get_point(i) for i in range(st.num_vertices())]) # We want to plot the alpha-complex with alpha=0.1. diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx index ea128743..5d181391 100644 --- a/src/python/gudhi/alpha_complex.pyx +++ b/src/python/gudhi/alpha_complex.pyx @@ -16,7 +16,9 @@ from libcpp.utility cimport pair from libcpp.string cimport string from libcpp cimport bool from libc.stdint cimport intptr_t +import errno import os +import warnings from gudhi.simplex_tree cimport * from gudhi.simplex_tree import SimplexTree @@ -28,7 +30,7 @@ __license__ = "GPL v3" cdef extern from "Alpha_complex_interface.h" namespace "Gudhi": cdef cppclass Alpha_complex_interface "Gudhi::alpha_complex::Alpha_complex_interface": - Alpha_complex_interface(vector[vector[double]] points, bool fast_version, bool exact_version) nogil except + + Alpha_complex_interface(vector[vector[double]] points, vector[double] weights, bool fast_version, bool exact_version) nogil except + vector[double] get_point(int vertex) nogil except + void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool default_filtration_value) nogil except + @@ -55,39 +57,56 @@ cdef class AlphaComplex: cdef Alpha_complex_interface * this_ptr # Fake constructor that does nothing but documenting the constructor - def __init__(self, points=None, off_file='', precision='safe'): + def __init__(self, points=[], off_file='', weights=[], precision='safe'): """AlphaComplex constructor. :param points: A list of points in d-Dimension. - :type points: list of list of double + :type points: Iterable[Iterable[float]] - Or - - :param off_file: An OFF file style name. + :param off_file: **[deprecated]** An `OFF file style <fileformats.html#off-file-format>`_ + name. + If an `off_file` is given with `points` as arguments, only points from the file are + taken into account. :type off_file: string - :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + :param weights: A list of weights. If set, the number of weights must correspond to the + number of points. + :type weights: Iterable[float] + + :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is + 'safe'. :type precision: string + + :raises FileNotFoundError: **[deprecated]** If `off_file` is set but not found. + :raises ValueError: In case of inconsistency between the number of points and weights. """ # The real cython constructor - def __cinit__(self, points = None, off_file = '', precision = 'safe'): - assert precision in ['fast', 'safe', 'exact'], "Alpha complex precision can only be 'fast', 'safe' or 'exact'" + def __cinit__(self, points = [], off_file = '', weights=[], precision = 'safe'): + assert precision in ['fast', 'safe', 'exact'], \ + "Alpha complex precision can only be 'fast', 'safe' or 'exact'" cdef bool fast = precision == 'fast' cdef bool exact = precision == 'exact' - cdef vector[vector[double]] pts if off_file: + warnings.warn("off_file is a deprecated parameter, please consider using gudhi.read_points_from_off_file", + DeprecationWarning) if os.path.isfile(off_file): points = read_points_from_off_file(off_file = off_file) else: - print("file " + off_file + " not found.") - if points is None: - # Empty Alpha construction - points=[] + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), off_file) + + # weights are set but is inconsistent with the number of points + if len(weights) != 0 and len(weights) != len(points): + raise ValueError("Inconsistency between the number of points and weights") + + # need to copy the points to use them without the gil + cdef vector[vector[double]] pts + cdef vector[double] wgts pts = points + wgts = weights with nogil: - self.this_ptr = new Alpha_complex_interface(pts, fast, exact) + self.this_ptr = new Alpha_complex_interface(pts, wgts, fast, exact) def __dealloc__(self): if self.this_ptr != NULL: diff --git a/src/python/gudhi/alpha_complex_3d.pyx b/src/python/gudhi/alpha_complex_3d.pyx new file mode 100644 index 00000000..3959004a --- /dev/null +++ b/src/python/gudhi/alpha_complex_3d.pyx @@ -0,0 +1,129 @@ +# 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 + +from __future__ import print_function +from cython cimport numeric +from libcpp.vector cimport vector +from libcpp.utility cimport pair +from libcpp.string cimport string +from libcpp cimport bool +from libc.stdint cimport intptr_t +import errno +import os +import warnings + +from gudhi.simplex_tree cimport * +from gudhi.simplex_tree import SimplexTree +from gudhi import read_points_from_off_file + +__author__ = "Vincent Rouvreau" +__copyright__ = "Copyright (C) 2021 Inria" +__license__ = "GPL v3" + +cdef extern from "Alpha_complex_interface_3d.h" namespace "Gudhi": + cdef cppclass Alpha_complex_interface_3d "Gudhi::alpha_complex::Alpha_complex_interface_3d": + Alpha_complex_interface_3d(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) nogil except + + +# AlphaComplex3D python interface +cdef class AlphaComplex3D: + """AlphaComplex3D is a simplicial complex constructed from the finite cells + of a Delaunay Triangulation. + + The filtration value of each simplex is computed as the square of the + circumradius of the simplex if the circumsphere is empty (the simplex is + then said to be Gabriel), and as the minimum of the filtration values of + the codimension 1 cofaces that make it not Gabriel otherwise. + + All simplices that have a filtration value strictly greater than a given + alpha squared value are not inserted into the complex. + + .. note:: + + When AlphaComplex3D is constructed with an infinite value of alpha, the + complex is a Delaunay complex. + + """ + + cdef Alpha_complex_interface_3d * this_ptr + + # Fake constructor that does nothing but documenting the constructor + def __init__(self, points=[], weights=[], precision='safe'): + """AlphaComplex3D constructor. + + :param points: A list of points in d-Dimension. + :type points: Iterable[Iterable[float]] + + :param weights: A list of weights. If set, the number of weights must correspond to the + number of points. + :type weights: Iterable[float] + + :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is + 'safe'. + :type precision: string + + :raises ValueError: In case of inconsistency between the number of points and weights. + """ + + # The real cython constructor + def __cinit__(self, points = [], weights=[], precision = 'safe'): + assert precision in ['fast', 'safe', 'exact'], \ + "Alpha complex precision can only be 'fast', 'safe' or 'exact'" + cdef bool fast = precision == 'fast' + cdef bool exact = precision == 'exact' + + # weights are set but is inconsistent with the number of points + if len(weights) != 0 and len(weights) != len(points): + raise ValueError("Inconsistency between the number of points and weights") + + # need to copy the points to use them without the gil + cdef vector[vector[double]] pts + cdef vector[double] wgts + pts = points + wgts = weights + with nogil: + self.this_ptr = new Alpha_complex_interface_3d(pts, wgts, fast, exact) + + def __dealloc__(self): + if self.this_ptr != NULL: + del self.this_ptr + + def __is_defined(self): + """Returns true if AlphaComplex3D pointer is not NULL. + """ + return self.this_ptr != NULL + + def get_point(self, vertex): + """This function returns the point corresponding to a given vertex from the :class:`~gudhi.SimplexTree`. + + :param vertex: The vertex. + :type vertex: int + :rtype: list of float + :returns: the point. + """ + return self.this_ptr.get_point(vertex) + + def create_simplex_tree(self, max_alpha_square = float('inf')): + """ + :param max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to + infinity, and there is very little point using anything else since it does not save time. + :type max_alpha_square: float + :returns: A simplex tree created from the Delaunay Triangulation. + :rtype: SimplexTree + """ + stree = SimplexTree() + cdef double mas = max_alpha_square + cdef intptr_t stree_int_ptr=stree.thisptr + with nogil: + self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr, + mas) + return stree diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index 3405fdd6..7d45af7c 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -31,15 +31,38 @@ namespace Gudhi { namespace alpha_complex { -template <typename CgalPointType> -std::vector<double> pt_cgal_to_cython(CgalPointType const& point) { - std::vector<double> vd; - vd.reserve(point.dimension()); - for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) - vd.push_back(CGAL::to_double(*coord)); - return vd; -} +// template Functor that transforms a CGAL point to a vector of double as expected by cython +template<typename CgalPointType, bool Weighted> +struct Point_cgal_to_cython; + +// Specialized Unweighted Functor +template<typename CgalPointType> +struct Point_cgal_to_cython<CgalPointType, false> { + std::vector<double> operator()(CgalPointType const& point) const + { + std::vector<double> vd; + vd.reserve(point.dimension()); + for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) + vd.push_back(CGAL::to_double(*coord)); + return vd; + } +}; +// Specialized Weighted Functor +template<typename CgalPointType> +struct Point_cgal_to_cython<CgalPointType, true> { + std::vector<double> operator()(CgalPointType const& weighted_point) const + { + auto point = weighted_point.point(); + std::vector<double> vd; + vd.reserve(point.dimension()); + for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) + vd.push_back(CGAL::to_double(*coord)); + return vd; + } +}; + +// Function that transforms a cython point (aka. a vector of double) to a CGAL point template <typename CgalPointType> static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) { return CgalPointType(vec.size(), vec.begin(), vec.end()); @@ -55,20 +78,29 @@ class Abstract_alpha_complex { virtual ~Abstract_alpha_complex() = default; }; -class Exact_Alphacomplex_dD final : public Abstract_alpha_complex { +template <bool Weighted = false> +class Exact_alpha_complex_dD final : public Abstract_alpha_complex { private: using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; - using Point = typename Kernel::Point_d; + using Bare_point = typename Kernel::Point_d; + using Point = std::conditional_t<Weighted, typename Kernel::Weighted_point_d, + typename Kernel::Point_d>; public: - Exact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version) + Exact_alpha_complex_dD(const std::vector<std::vector<double>>& points, bool exact_version) : exact_version_(exact_version), - alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>)) { + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>)) { + } + + Exact_alpha_complex_dD(const std::vector<std::vector<double>>& points, + const std::vector<double>& weights, bool exact_version) + : exact_version_(exact_version), + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>), weights) { } virtual std::vector<double> get_point(int vh) override { - Point const& point = alpha_complex_.get_point(vh); - return pt_cgal_to_cython(point); + // Can be a Weighted or a Bare point in function of Weighted + return Point_cgal_to_cython<Point, Weighted>()(alpha_complex_.get_point(vh)); } virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, @@ -78,23 +110,32 @@ class Exact_Alphacomplex_dD final : public Abstract_alpha_complex { private: bool exact_version_; - Alpha_complex<Kernel> alpha_complex_; + Alpha_complex<Kernel, Weighted> alpha_complex_; }; -class Inexact_Alphacomplex_dD final : public Abstract_alpha_complex { +template <bool Weighted = false> +class Inexact_alpha_complex_dD final : public Abstract_alpha_complex { private: using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; - using Point = typename Kernel::Point_d; + using Bare_point = typename Kernel::Point_d; + using Point = std::conditional_t<Weighted, typename Kernel::Weighted_point_d, + typename Kernel::Point_d>; public: - Inexact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version) + Inexact_alpha_complex_dD(const std::vector<std::vector<double>>& points, bool exact_version) : exact_version_(exact_version), - alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>)) { + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>)) { + } + + Inexact_alpha_complex_dD(const std::vector<std::vector<double>>& points, + const std::vector<double>& weights, bool exact_version) + : exact_version_(exact_version), + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Bare_point>), weights) { } virtual std::vector<double> get_point(int vh) override { - Point const& point = alpha_complex_.get_point(vh); - return pt_cgal_to_cython(point); + // Can be a Weighted or a Bare point in function of Weighted + return Point_cgal_to_cython<Point, Weighted>()(alpha_complex_.get_point(vh)); } virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) override { @@ -103,35 +144,41 @@ class Inexact_Alphacomplex_dD final : public Abstract_alpha_complex { private: bool exact_version_; - Alpha_complex<Kernel> alpha_complex_; + Alpha_complex<Kernel, Weighted> alpha_complex_; }; -template <complexity Complexity> -class Alphacomplex_3D final : public Abstract_alpha_complex { +template <complexity Complexity, bool Weighted = false> +class Alpha_complex_3D final : public Abstract_alpha_complex { private: - using Point = typename Alpha_complex_3d<Complexity, false, false>::Bare_point_3; + using Bare_point = typename Alpha_complex_3d<Complexity, Weighted, false>::Bare_point_3; + using Point = typename Alpha_complex_3d<Complexity, Weighted, false>::Point_3; - static Point pt_cython_to_cgal_3(std::vector<double> const& vec) { - return Point(vec[0], vec[1], vec[2]); + static Bare_point pt_cython_to_cgal_3(std::vector<double> const& vec) { + return Bare_point(vec[0], vec[1], vec[2]); } public: - Alphacomplex_3D(const std::vector<std::vector<double>>& points) + Alpha_complex_3D(const std::vector<std::vector<double>>& points) : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3)) { } + Alpha_complex_3D(const std::vector<std::vector<double>>& points, const std::vector<double>& weights) + : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3), weights) { + } + virtual std::vector<double> get_point(int vh) override { - Point const& point = alpha_complex_.get_point(vh); - return pt_cgal_to_cython(point); + // Can be a Weighted or a Bare point in function of Weighted + return Point_cgal_to_cython<Point, Weighted>()(alpha_complex_.get_point(vh)); } virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) override { - return alpha_complex_.create_complex(*simplex_tree, max_alpha_square); + alpha_complex_.create_complex(*simplex_tree, max_alpha_square); + return true; } private: - Alpha_complex_3d<Complexity, false, false> alpha_complex_; + Alpha_complex_3d<Complexity, Weighted, false> alpha_complex_; }; diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 23be194d..ed243f19 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -27,10 +27,24 @@ namespace alpha_complex { class Alpha_complex_interface { public: - Alpha_complex_interface(const std::vector<std::vector<double>>& points, bool fast_version, bool exact_version) - : points_(points), - fast_version_(fast_version), - exact_version_(exact_version) { + Alpha_complex_interface(const std::vector<std::vector<double>>& points, + const std::vector<double>& weights, + bool fast_version, bool exact_version) + : empty_point_set_(points.size() == 0) { + const bool weighted = (weights.size() > 0); + if (fast_version) { + if (weighted) { + alpha_ptr_ = std::make_unique<Inexact_alpha_complex_dD<true>>(points, weights, exact_version); + } else { + alpha_ptr_ = std::make_unique<Inexact_alpha_complex_dD<false>>(points, exact_version); + } + } else { + if (weighted) { + alpha_ptr_ = std::make_unique<Exact_alpha_complex_dD<true>>(points, weights, exact_version); + } else { + alpha_ptr_ = std::make_unique<Exact_alpha_complex_dD<false>>(points, exact_version); + } + } } std::vector<double> get_point(int vh) { @@ -39,38 +53,14 @@ class Alpha_complex_interface { void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) { - if (points_.size() > 0) { - std::size_t dimension = points_[0].size(); - if (dimension == 3 && !default_filtration_value) { - if (fast_version_) - alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::FAST>>(points_); - else if (exact_version_) - alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::EXACT>>(points_); - else - alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::SAFE>>(points_); - if (!alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value)) { - // create_simplex_tree will fail if all points are on a plane - Retry with dD by setting dimension to 2 - dimension--; - alpha_ptr_.reset(); - } - } - // Not ** else ** because we have to take into account if 3d fails - if (dimension != 3 || default_filtration_value) { - if (fast_version_) { - alpha_ptr_ = std::make_unique<Inexact_Alphacomplex_dD>(points_, exact_version_); - } else { - alpha_ptr_ = std::make_unique<Exact_Alphacomplex_dD>(points_, exact_version_); - } - alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); - } - } + // Nothing to be done in case of an empty point set + if (!empty_point_set_) + alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); } private: std::unique_ptr<Abstract_alpha_complex> alpha_ptr_; - std::vector<std::vector<double>> points_; - bool fast_version_; - bool exact_version_; + bool empty_point_set_; }; } // namespace alpha_complex diff --git a/src/python/include/Alpha_complex_interface_3d.h b/src/python/include/Alpha_complex_interface_3d.h new file mode 100644 index 00000000..bb66b8e1 --- /dev/null +++ b/src/python/include/Alpha_complex_interface_3d.h @@ -0,0 +1,71 @@ +/* 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 + */ + +#ifndef INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ +#define INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ + +#include "Alpha_complex_factory.h" +#include <gudhi/Alpha_complex_options.h> + +#include "Simplex_tree_interface.h" + +#include <iostream> +#include <vector> +#include <string> +#include <memory> // for std::unique_ptr + +namespace Gudhi { + +namespace alpha_complex { + +class Alpha_complex_interface_3d { + public: + Alpha_complex_interface_3d(const std::vector<std::vector<double>>& points, + const std::vector<double>& weights, + bool fast_version, bool exact_version) + : empty_point_set_(points.size() == 0) { + const bool weighted = (weights.size() > 0); + if (fast_version) + if (weighted) + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::FAST, true>>(points, weights); + else + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::FAST>>(points); + else if (exact_version) + if (weighted) + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::EXACT, true>>(points, weights); + else + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::EXACT>>(points); + else + if (weighted) + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::SAFE, true>>(points, weights); + else + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::SAFE>>(points); + } + + std::vector<double> get_point(int vh) { + return alpha_ptr_->get_point(vh); + } + + void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square) { + // Nothing to be done in case of an empty point set + if (!empty_point_set_) + alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, false); + } + + private: + std::unique_ptr<Abstract_alpha_complex> alpha_ptr_; + bool empty_point_set_; +}; + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py index 814f8289..e0f2b5df 100755 --- a/src/python/test/test_alpha_complex.py +++ b/src/python/test/test_alpha_complex.py @@ -12,6 +12,8 @@ import gudhi as gd import math import numpy as np import pytest +import warnings + try: # python3 from itertools import zip_longest @@ -25,12 +27,17 @@ __license__ = "MIT" def _empty_alpha(precision): + alpha_complex = gd.AlphaComplex(precision = precision) + assert alpha_complex.__is_defined() == True + +def _one_2d_point_alpha(precision): alpha_complex = gd.AlphaComplex(points=[[0, 0]], precision = precision) assert alpha_complex.__is_defined() == True def test_empty_alpha(): for precision in ['fast', 'safe', 'exact']: _empty_alpha(precision) + _one_2d_point_alpha(precision) def _infinite_alpha(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] @@ -198,7 +205,13 @@ def test_delaunay_complex(): _delaunay_complex(precision) def _3d_points_on_a_plane(precision, default_filtration_value): - alpha = gd.AlphaComplex(off_file='alphacomplexdoc.off', precision = precision) + alpha = gd.AlphaComplex(points = [[1.0, 1.0 , 0.0], + [7.0, 0.0 , 0.0], + [4.0, 6.0 , 0.0], + [9.0, 6.0 , 0.0], + [0.0, 14.0, 0.0], + [2.0, 19.0, 0.0], + [9.0, 17.0, 0.0]], precision = precision) simplex_tree = alpha.create_simplex_tree(default_filtration_value = default_filtration_value) assert simplex_tree.dimension() == 2 @@ -206,28 +219,16 @@ def _3d_points_on_a_plane(precision, default_filtration_value): assert simplex_tree.num_simplices() == 25 def test_3d_points_on_a_plane(): - off_file = open("alphacomplexdoc.off", "w") - off_file.write("OFF \n" \ - "7 0 0 \n" \ - "1.0 1.0 0.0\n" \ - "7.0 0.0 0.0\n" \ - "4.0 6.0 0.0\n" \ - "9.0 6.0 0.0\n" \ - "0.0 14.0 0.0\n" \ - "2.0 19.0 0.0\n" \ - "9.0 17.0 0.0\n" ) - off_file.close() - for default_filtration_value in [True, False]: for precision in ['fast', 'safe', 'exact']: _3d_points_on_a_plane(precision, default_filtration_value) def _3d_tetrahedrons(precision): points = 10*np.random.rand(10, 3) - alpha = gd.AlphaComplex(points=points, precision = precision) + alpha = gd.AlphaComplex(points = points, precision = precision) st_alpha = alpha.create_simplex_tree(default_filtration_value = False) # New AlphaComplex for get_point to work - delaunay = gd.AlphaComplex(points=points, precision = precision) + delaunay = gd.AlphaComplex(points = points, precision = precision) st_delaunay = delaunay.create_simplex_tree(default_filtration_value = True) delaunay_tetra = [] @@ -256,3 +257,59 @@ def _3d_tetrahedrons(precision): def test_3d_tetrahedrons(): for precision in ['fast', 'safe', 'exact']: _3d_tetrahedrons(precision) + +def test_off_file_deprecation_warning(): + off_file = open("alphacomplexdoc.off", "w") + off_file.write("OFF \n" \ + "7 0 0 \n" \ + "1.0 1.0 0.0\n" \ + "7.0 0.0 0.0\n" \ + "4.0 6.0 0.0\n" \ + "9.0 6.0 0.0\n" \ + "0.0 14.0 0.0\n" \ + "2.0 19.0 0.0\n" \ + "9.0 17.0 0.0\n" ) + off_file.close() + + with pytest.warns(DeprecationWarning): + alpha = gd.AlphaComplex(off_file="alphacomplexdoc.off") + +def test_non_existing_off_file(): + with pytest.raises(FileNotFoundError): + alpha = gd.AlphaComplex(off_file="pouetpouettralala.toubiloubabdou") + +def test_inconsistency_points_and_weights(): + points = [[1.0, 1.0 , 0.0], + [7.0, 0.0 , 0.0], + [4.0, 6.0 , 0.0], + [9.0, 6.0 , 0.0], + [0.0, 14.0, 0.0], + [2.0, 19.0, 0.0], + [9.0, 17.0, 0.0]] + with pytest.raises(ValueError): + # 7 points, 8 weights, on purpose + alpha = gd.AlphaComplex(points = points, + weights = [1., 2., 3., 4., 5., 6., 7., 8.]) + + with pytest.raises(ValueError): + # 7 points, 6 weights, on purpose + alpha = gd.AlphaComplex(points = points, + weights = [1., 2., 3., 4., 5., 6.]) + +def _doc_example(precision): + stree_from_values = gd.AlphaComplex(points=[[ 1., -1., -1.], + [-1., 1., -1.], + [-1., -1., 1.], + [ 1., 1., 1.], + [ 2., 2., 2.]], + weights = [4., 4., 4., 4., 1.], + precision = precision).create_simplex_tree() + + assert stree_from_values.filtration([0, 1, 2, 3]) == pytest.approx(-1.) + assert stree_from_values.filtration([0, 1, 3, 4]) == pytest.approx(95.) + assert stree_from_values.filtration([0, 2, 3, 4]) == pytest.approx(95.) + assert stree_from_values.filtration([1, 2, 3, 4]) == pytest.approx(95.) + +def test_doc_example(): + for precision in ['fast', 'safe', 'exact']: + _doc_example(precision) diff --git a/src/python/test/test_reader_utils.py b/src/python/test/test_reader_utils.py index 90da6651..4fc7c00f 100755 --- a/src/python/test/test_reader_utils.py +++ b/src/python/test/test_reader_utils.py @@ -8,8 +8,9 @@ - YYYY/MM Author: Description of the modification """ -import gudhi +import gudhi as gd import numpy as np +from pytest import raises __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2017 Inria" @@ -18,7 +19,7 @@ __license__ = "MIT" def test_non_existing_csv_file(): # Try to open a non existing file - matrix = gudhi.read_lower_triangular_matrix_from_csv_file( + matrix = gd.read_lower_triangular_matrix_from_csv_file( csv_file="pouetpouettralala.toubiloubabdou" ) assert matrix == [] @@ -29,7 +30,7 @@ def test_full_square_distance_matrix_csv_file(): test_file = open("full_square_distance_matrix.csv", "w") test_file.write("0;1;2;3;\n1;0;4;5;\n2;4;0;6;\n3;5;6;0;") test_file.close() - matrix = gudhi.read_lower_triangular_matrix_from_csv_file( + matrix = gd.read_lower_triangular_matrix_from_csv_file( csv_file="full_square_distance_matrix.csv" ) assert matrix == [[], [1.0], [2.0, 4.0], [3.0, 5.0, 6.0]] @@ -40,7 +41,7 @@ def test_lower_triangular_distance_matrix_csv_file(): test_file = open("lower_triangular_distance_matrix.csv", "w") test_file.write("\n1,\n2,3,\n4,5,6,\n7,8,9,10,") test_file.close() - matrix = gudhi.read_lower_triangular_matrix_from_csv_file( + matrix = gd.read_lower_triangular_matrix_from_csv_file( csv_file="lower_triangular_distance_matrix.csv", separator="," ) assert matrix == [[], [1.0], [2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0, 10.0]] @@ -48,11 +49,11 @@ def test_lower_triangular_distance_matrix_csv_file(): def test_non_existing_persistence_file(): # Try to open a non existing file - persistence = gudhi.read_persistence_intervals_grouped_by_dimension( + persistence = gd.read_persistence_intervals_grouped_by_dimension( persistence_file="pouetpouettralala.toubiloubabdou" ) assert persistence == [] - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="pouetpouettralala.toubiloubabdou", only_this_dim=1 ) np.testing.assert_array_equal(persistence, []) @@ -65,21 +66,21 @@ def test_read_persistence_intervals_without_dimension(): "# Simple persistence diagram without dimension\n2.7 3.7\n9.6 14.\n34.2 34.974\n3. inf" ) test_file.close() - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_without_dimension.pers" ) np.testing.assert_array_equal( persistence, [(2.7, 3.7), (9.6, 14.0), (34.2, 34.974), (3.0, float("Inf"))] ) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_without_dimension.pers", only_this_dim=0 ) np.testing.assert_array_equal(persistence, []) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_without_dimension.pers", only_this_dim=1 ) np.testing.assert_array_equal(persistence, []) - persistence = gudhi.read_persistence_intervals_grouped_by_dimension( + persistence = gd.read_persistence_intervals_grouped_by_dimension( persistence_file="persistence_intervals_without_dimension.pers" ) assert persistence == { @@ -94,29 +95,29 @@ def test_read_persistence_intervals_with_dimension(): "# Simple persistence diagram with dimension\n0 2.7 3.7\n1 9.6 14.\n3 34.2 34.974\n1 3. inf" ) test_file.close() - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_with_dimension.pers" ) np.testing.assert_array_equal( persistence, [(2.7, 3.7), (9.6, 14.0), (34.2, 34.974), (3.0, float("Inf"))] ) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=0 ) np.testing.assert_array_equal(persistence, [(2.7, 3.7)]) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=1 ) np.testing.assert_array_equal(persistence, [(9.6, 14.0), (3.0, float("Inf"))]) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=2 ) np.testing.assert_array_equal(persistence, []) - persistence = gudhi.read_persistence_intervals_in_dimension( + persistence = gd.read_persistence_intervals_in_dimension( persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=3 ) np.testing.assert_array_equal(persistence, [(34.2, 34.974)]) - persistence = gudhi.read_persistence_intervals_grouped_by_dimension( + persistence = gd.read_persistence_intervals_grouped_by_dimension( persistence_file="persistence_intervals_with_dimension.pers" ) assert persistence == { |