diff options
author | Vincent Rouvreau <10407034+VincentRouvreau@users.noreply.github.com> | 2020-06-06 14:44:14 +0200 |
---|---|---|
committer | GitHub <noreply@github.com> | 2020-06-06 14:44:14 +0200 |
commit | 434cc66fbf832ca492e5f0710c8afc67d681637a (patch) | |
tree | 08cc98efa6a60aeaa59cf6407d7c4e722dbb212b /src/python | |
parent | d1df63c72af9a239b5dad44d0f8545d24a4a292b (diff) | |
parent | 19dc69e6edaac076da5138fe091be367a3652057 (diff) |
Merge pull request #326 from VincentRouvreau/alpha_complex_python_fast_and_exact
Alpha complex python enhancement
Diffstat (limited to 'src/python')
-rw-r--r-- | src/python/doc/alpha_complex_user.rst | 20 | ||||
-rw-r--r-- | src/python/gudhi/alpha_complex.pyx | 49 | ||||
-rw-r--r-- | src/python/include/Alpha_complex_interface.h | 69 | ||||
-rwxr-xr-x | src/python/test/test_alpha_complex.py | 75 |
4 files changed, 162 insertions, 51 deletions
diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst index e31194a7..097470c1 100644 --- a/src/python/doc/alpha_complex_user.rst +++ b/src/python/doc/alpha_complex_user.rst @@ -16,8 +16,24 @@ Definition Remarks ^^^^^^^ -When an :math:`\alpha`-complex is constructed with an infinite value of :math:`\alpha^2`, -the complex is a Delaunay complex (with special filtration values). +* When an :math:`\alpha`-complex is constructed with an infinite value of :math:`\alpha^2`, the complex is a Delaunay + complex (with special filtration values). The Delaunay complex without filtration values is also available by + passing :code:`default_filtration_value = True` to :func:`~gudhi.AlphaComplex.create_simplex_tree`. +* For people only interested in the topology of the Alpha complex (for instance persistence), Alpha complex is + equivalent to the `Čech complex <https://gudhi.inria.fr/doc/latest/group__cech__complex.html>`_ and much smaller if + you do not bound the radii. `Čech complex <https://gudhi.inria.fr/doc/latest/group__cech__complex.html>`_ can still + make sense in higher dimension precisely because you can bound the radii. +* Using the default :code:`precision = 'safe'` makes the construction safe. + 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. + 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 + the computation of filtration values can exceptionally be arbitrarily bad. In all cases, we still guarantee that the + output is a valid filtration (faces have a filtration value no larger than their cofaces). +* For performances reasons, it is advised to use Alpha_complex with `CGAL <installation.html#cgal>`_ :math:`\geq` 5.0.0. For performances reasons, it is advised to use CGAL :math:`\geq` 5.0.0. diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx index d75e374a..a356384d 100644 --- a/src/python/gudhi/alpha_complex.pyx +++ b/src/python/gudhi/alpha_complex.pyx @@ -27,11 +27,11 @@ __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) nogil except + + Alpha_complex_interface(vector[vector[double]] points, bool fast_version) nogil except + # bool from_file is a workaround for cython to find the correct signature - Alpha_complex_interface(string off_file, bool from_file) nogil except + + Alpha_complex_interface(string off_file, bool fast_version, bool from_file) 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 + + void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool exact_version, bool default_filtration_value) nogil except + # AlphaComplex python interface cdef class AlphaComplex: @@ -53,10 +53,11 @@ cdef class AlphaComplex: """ - cdef Alpha_complex_interface * thisptr + cdef Alpha_complex_interface * this_ptr + cdef bool exact # Fake constructor that does nothing but documenting the constructor - def __init__(self, points=None, off_file=''): + def __init__(self, points=None, off_file='', precision='safe'): """AlphaComplex constructor. :param points: A list of points in d-Dimension. @@ -66,15 +67,21 @@ cdef class AlphaComplex: :param off_file: An OFF file style name. :type off_file: string + + :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + :type precision: string """ # The real cython constructor - def __cinit__(self, points = None, off_file = ''): + 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'" + cdef bool fast = precision == 'fast' + self.exact = precision == 'exact' + cdef vector[vector[double]] pts if off_file: if os.path.isfile(off_file): - self.thisptr = new Alpha_complex_interface( - off_file.encode('utf-8'), True) + self.this_ptr = new Alpha_complex_interface(off_file.encode('utf-8'), fast, True) else: print("file " + off_file + " not found.") else: @@ -83,17 +90,16 @@ cdef class AlphaComplex: points=[] pts = points with nogil: - self.thisptr = new Alpha_complex_interface(pts) - + self.this_ptr = new Alpha_complex_interface(pts, fast) def __dealloc__(self): - if self.thisptr != NULL: - del self.thisptr + if self.this_ptr != NULL: + del self.this_ptr def __is_defined(self): """Returns true if AlphaComplex pointer is not NULL. """ - return self.thisptr != NULL + return self.this_ptr != NULL def get_point(self, vertex): """This function returns the point corresponding to a given vertex. @@ -103,21 +109,24 @@ cdef class AlphaComplex: :rtype: list of float :returns: the point. """ - return self.thisptr.get_point(vertex) + return self.this_ptr.get_point(vertex) - def create_simplex_tree(self, max_alpha_square = float('inf')): + def create_simplex_tree(self, max_alpha_square = float('inf'), default_filtration_value = False): """ - :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. + :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 + :param default_filtration_value: Set this value to `True` if filtration values are not needed to be computed + (will be set to `NaN`). Default value is `False` (which means compute the filtration values). + :type default_filtration_value: bool :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 + cdef bool compute_filtration = default_filtration_value == True with nogil: - self.thisptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr, mas) + self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr, + mas, self.exact, compute_filtration) return stree diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 40de88f3..3ac5db1f 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -23,45 +23,74 @@ #include <iostream> #include <vector> #include <string> +#include <memory> // for std::unique_ptr namespace Gudhi { namespace alpha_complex { class Alpha_complex_interface { - using Dynamic_kernel = CGAL::Epeck_d< CGAL::Dynamic_dimension_tag >; - using Point_d = Dynamic_kernel::Point_d; + private: + using Exact_kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; + using Inexact_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; + using Point_exact_kernel = typename Exact_kernel::Point_d; + using Point_inexact_kernel = typename Inexact_kernel::Point_d; - public: - Alpha_complex_interface(const std::vector<std::vector<double>>& points) { - auto mkpt = [](std::vector<double> const& vec){ - return Point_d(vec.size(), vec.begin(), vec.end()); - }; - alpha_complex_ = new Alpha_complex<Dynamic_kernel>(boost::adaptors::transform(points, mkpt)); + template <typename CgalPointType> + std::vector<double> pt_cgal_to_cython(CgalPointType& point) { + std::vector<double> vd; + for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) + vd.push_back(CGAL::to_double(*coord)); + return vd; } - Alpha_complex_interface(const std::string& off_file_name, bool from_file = true) { - alpha_complex_ = new Alpha_complex<Dynamic_kernel>(off_file_name); + template <typename CgalPointType> + static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) { + return CgalPointType(vec.size(), vec.begin(), vec.end()); } - ~Alpha_complex_interface() { - delete alpha_complex_; + public: + Alpha_complex_interface(const std::vector<std::vector<double>>& points, bool fast_version) + : fast_version_(fast_version) { + if (fast_version_) { + ac_inexact_ptr_ = std::make_unique<Alpha_complex<Inexact_kernel>>( + boost::adaptors::transform(points, pt_cython_to_cgal<Point_inexact_kernel>)); + } else { + ac_exact_ptr_ = std::make_unique<Alpha_complex<Exact_kernel>>( + boost::adaptors::transform(points, pt_cython_to_cgal<Point_exact_kernel>)); + } + } + + Alpha_complex_interface(const std::string& off_file_name, bool fast_version, bool from_file = true) + : fast_version_(fast_version) { + if (fast_version_) + ac_inexact_ptr_ = std::make_unique<Alpha_complex<Inexact_kernel>>(off_file_name); + else + ac_exact_ptr_ = std::make_unique<Alpha_complex<Exact_kernel>>(off_file_name); } std::vector<double> get_point(int vh) { - std::vector<double> vd; - Point_d const& ph = alpha_complex_->get_point(vh); - for (auto coord = ph.cartesian_begin(); coord != ph.cartesian_end(); coord++) - vd.push_back(CGAL::to_double(*coord)); - return vd; + if (fast_version_) { + Point_inexact_kernel const& point = ac_inexact_ptr_->get_point(vh); + return pt_cgal_to_cython(point); + } else { + Point_exact_kernel const& point = ac_exact_ptr_->get_point(vh); + return pt_cgal_to_cython(point); + } } - void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square) { - alpha_complex_->create_complex(*simplex_tree, max_alpha_square); + void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool exact_version, + bool default_filtration_value) { + if (fast_version_) + ac_inexact_ptr_->create_complex(*simplex_tree, max_alpha_square, exact_version, default_filtration_value); + else + ac_exact_ptr_->create_complex(*simplex_tree, max_alpha_square, exact_version, default_filtration_value); } private: - Alpha_complex<Dynamic_kernel>* alpha_complex_; + bool fast_version_; + std::unique_ptr<Alpha_complex<Exact_kernel>> ac_exact_ptr_; + std::unique_ptr<Alpha_complex<Inexact_kernel>> ac_inexact_ptr_; }; } // namespace alpha_complex diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py index 77121302..943ad2c4 100755 --- a/src/python/test/test_alpha_complex.py +++ b/src/python/test/test_alpha_complex.py @@ -24,14 +24,18 @@ __copyright__ = "Copyright (C) 2016 Inria" __license__ = "MIT" -def test_empty_alpha(): - alpha_complex = AlphaComplex(points=[[0, 0]]) +def _empty_alpha(precision): + alpha_complex = AlphaComplex(points=[[0, 0]], precision = precision) assert alpha_complex.__is_defined() == True +def test_empty_alpha(): + _empty_alpha('fast') + _empty_alpha('safe') + _empty_alpha('exact') -def test_infinite_alpha(): +def _infinite_alpha(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - alpha_complex = AlphaComplex(points=point_list) + alpha_complex = AlphaComplex(points=point_list, precision = precision) assert alpha_complex.__is_defined() == True simplex_tree = alpha_complex.create_simplex_tree() @@ -79,10 +83,14 @@ def test_infinite_alpha(): else: assert False +def test_infinite_alpha(): + _infinite_alpha('fast') + _infinite_alpha('safe') + _infinite_alpha('exact') -def test_filtered_alpha(): +def _filtered_alpha(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - filtered_alpha = AlphaComplex(points=point_list) + filtered_alpha = AlphaComplex(points=point_list, precision = precision) simplex_tree = filtered_alpha.create_simplex_tree(max_alpha_square=0.25) @@ -119,7 +127,12 @@ def test_filtered_alpha(): assert simplex_tree.get_star([0]) == [([0], 0.0), ([0, 1], 0.25), ([0, 2], 0.25)] assert simplex_tree.get_cofaces([0], 1) == [([0, 1], 0.25), ([0, 2], 0.25)] -def test_safe_alpha_persistence_comparison(): +def test_filtered_alpha(): + _filtered_alpha('fast') + _filtered_alpha('safe') + _filtered_alpha('exact') + +def _safe_alpha_persistence_comparison(precision): #generate periodic signal time = np.arange(0, 10, 1) signal = [math.sin(x) for x in time] @@ -131,10 +144,10 @@ def test_safe_alpha_persistence_comparison(): embedding2 = [[signal[i], delayed[i]] for i in range(len(time))] #build alpha complex and simplex tree - alpha_complex1 = AlphaComplex(points=embedding1) + alpha_complex1 = AlphaComplex(points=embedding1, precision = precision) simplex_tree1 = alpha_complex1.create_simplex_tree() - alpha_complex2 = AlphaComplex(points=embedding2) + alpha_complex2 = AlphaComplex(points=embedding2, precision = precision) simplex_tree2 = alpha_complex2.create_simplex_tree() diag1 = simplex_tree1.persistence() @@ -143,3 +156,47 @@ def test_safe_alpha_persistence_comparison(): for (first_p, second_p) in zip_longest(diag1, diag2): assert first_p[0] == pytest.approx(second_p[0]) assert first_p[1] == pytest.approx(second_p[1]) + + +def test_safe_alpha_persistence_comparison(): + # Won't work for 'fast' version + _safe_alpha_persistence_comparison('safe') + _safe_alpha_persistence_comparison('exact') + +def _delaunay_complex(precision): + point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] + filtered_alpha = AlphaComplex(points=point_list, precision = precision) + + simplex_tree = filtered_alpha.create_simplex_tree(default_filtration_value = True) + + assert simplex_tree.num_simplices() == 11 + assert simplex_tree.num_vertices() == 4 + + assert point_list[0] == filtered_alpha.get_point(0) + assert point_list[1] == filtered_alpha.get_point(1) + assert point_list[2] == filtered_alpha.get_point(2) + assert point_list[3] == filtered_alpha.get_point(3) + try: + filtered_alpha.get_point(4) == [] + except IndexError: + pass + else: + assert False + try: + filtered_alpha.get_point(125) == [] + except IndexError: + pass + else: + assert False + + for filtered_value in simplex_tree.get_filtration(): + assert math.isnan(filtered_value[1]) + for filtered_value in simplex_tree.get_star([0]): + assert math.isnan(filtered_value[1]) + for filtered_value in simplex_tree.get_cofaces([0], 1): + assert math.isnan(filtered_value[1]) + +def test_delaunay_complex(): + _delaunay_complex('fast') + _delaunay_complex('safe') + _delaunay_complex('exact') |