diff options
Diffstat (limited to 'src/python')
-rw-r--r-- | src/python/doc/alpha_complex_user.rst | 47 | ||||
-rw-r--r-- | src/python/gudhi/alpha_complex.pyx | 27 | ||||
-rw-r--r-- | src/python/include/Alpha_complex_factory.h | 139 | ||||
-rw-r--r-- | src/python/include/Alpha_complex_interface.h | 92 | ||||
-rwxr-xr-x | src/python/test/test_alpha_complex.py | 83 |
5 files changed, 261 insertions, 127 deletions
diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst index 097470c1..bcc2fc51 100644 --- a/src/python/doc/alpha_complex_user.rst +++ b/src/python/doc/alpha_complex_user.rst @@ -178,49 +178,22 @@ In the following example, a threshold of :math:`\alpha^2 = 32.0` is used. Example from OFF file ^^^^^^^^^^^^^^^^^^^^^ -This example builds the Delaunay triangulation from the points given by an OFF file, and initializes the alpha complex -with it. +This example builds the alpha complex from 300 random points on a 2-torus. +Then, it computes the persistence diagram and diplays it: -Then, it is asked to display information about the alpha complex: - -.. testcode:: +.. plot:: + :include-source: + import matplotlib.pyplot as plt import gudhi alpha_complex = gudhi.AlphaComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/alphacomplexdoc.off') - simplex_tree = alpha_complex.create_simplex_tree(max_alpha_square=32.0) + '/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) - fmt = '%s -> %.2f' - for filtered_value in simplex_tree.get_filtration(): - print(fmt % tuple(filtered_value)) - -the program output is: - -.. testoutput:: - - Alpha complex is of dimension 2 - 20 simplices - 7 vertices. - [0] -> 0.00 - [1] -> 0.00 - [2] -> 0.00 - [3] -> 0.00 - [4] -> 0.00 - [5] -> 0.00 - [6] -> 0.00 - [2, 3] -> 6.25 - [4, 5] -> 7.25 - [0, 2] -> 8.50 - [0, 1] -> 9.25 - [1, 3] -> 10.00 - [1, 2] -> 11.25 - [1, 2, 3] -> 12.50 - [0, 1, 2] -> 13.00 - [5, 6] -> 13.25 - [2, 4] -> 20.00 - [4, 6] -> 22.74 - [4, 5, 6] -> 22.74 - [3, 6] -> 30.25 - + diag = simplex_tree.persistence() + gudhi.plot_persistence_diagram(diag) + plt.show() diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx index a356384d..ac44c61f 100644 --- a/src/python/gudhi/alpha_complex.pyx +++ b/src/python/gudhi/alpha_complex.pyx @@ -20,6 +20,7 @@ import os 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) 2016 Inria" @@ -27,11 +28,9 @@ __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) nogil except + - # bool from_file is a workaround for cython to find the correct signature - Alpha_complex_interface(string off_file, bool fast_version, bool from_file) nogil except + + Alpha_complex_interface(vector[vector[double]] points, 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 exact_version, bool default_filtration_value) nogil except + + void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool default_filtration_value) nogil except + # AlphaComplex python interface cdef class AlphaComplex: @@ -54,7 +53,6 @@ cdef class AlphaComplex: """ 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='', precision='safe'): @@ -76,21 +74,20 @@ cdef class AlphaComplex: 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 bool exact = precision == 'exact' cdef vector[vector[double]] pts if off_file: if os.path.isfile(off_file): - self.this_ptr = new Alpha_complex_interface(off_file.encode('utf-8'), fast, True) + points = read_points_from_off_file(off_file = off_file) else: print("file " + off_file + " not found.") - else: - if points is None: - # Empty Alpha construction - points=[] - pts = points - with nogil: - self.this_ptr = new Alpha_complex_interface(pts, fast) + if points is None: + # Empty Alpha construction + points=[] + pts = points + with nogil: + self.this_ptr = new Alpha_complex_interface(pts, fast, exact) def __dealloc__(self): if self.this_ptr != NULL: @@ -128,5 +125,5 @@ cdef class AlphaComplex: cdef bool compute_filtration = default_filtration_value == True with nogil: self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr, - mas, self.exact, compute_filtration) + mas, compute_filtration) return stree diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h new file mode 100644 index 00000000..7d98f454 --- /dev/null +++ b/src/python/include/Alpha_complex_factory.h @@ -0,0 +1,139 @@ +/* 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) 2020 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef INCLUDE_ALPHA_COMPLEX_FACTORY_H_ +#define INCLUDE_ALPHA_COMPLEX_FACTORY_H_ + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Alpha_complex.h> +#include <gudhi/Alpha_complex_3d.h> +#include <gudhi/Alpha_complex_options.h> +#include <CGAL/Epeck_d.h> +#include <CGAL/Epick_d.h> + +#include <boost/range/adaptor/transformed.hpp> + +#include "Simplex_tree_interface.h" + +#include <iostream> +#include <vector> +#include <string> +#include <memory> // for std::unique_ptr + +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 <typename CgalPointType> +static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) { + return CgalPointType(vec.size(), vec.begin(), vec.end()); +} + +class Abstract_alpha_complex { + public: + virtual std::vector<double> get_point(int vh) = 0; + virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, + bool default_filtration_value) = 0; +}; + +class Exact_Alphacomplex_dD : public Abstract_alpha_complex { + private: + using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; + using Point = typename Kernel::Point_d; + + public: + Exact_Alphacomplex_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>)) { + } + + std::vector<double> get_point(int vh) { + Point const& point = alpha_complex_.get_point(vh); + return pt_cgal_to_cython(point); + } + + bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, + bool default_filtration_value){ + return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); + } + + private: + bool exact_version_; + Alpha_complex<Kernel> alpha_complex_; +}; + +class Inexact_Alphacomplex_dD : public Abstract_alpha_complex { + private: + using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; + using Point = typename Kernel::Point_d; + + public: + Inexact_Alphacomplex_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>)) { + } + + std::vector<double> get_point(int vh) { + Point const& point = alpha_complex_.get_point(vh); + return pt_cgal_to_cython(point); + } + bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, + bool default_filtration_value){ + return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); + } + + private: + bool exact_version_; + Alpha_complex<Kernel> alpha_complex_; +}; + +template <complexity Complexity> +class Alphacomplex_3D : public Abstract_alpha_complex { + private: + using Point = typename Alpha_complex_3d<Complexity, false, false>::Bare_point_3; + + static Point pt_cython_to_cgal_3(std::vector<double> const& vec) { + return Point(vec[0], vec[1], vec[2]); + } + + public: + Alphacomplex_3D(const std::vector<std::vector<double>>& points) + : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3)) { + } + + std::vector<double> get_point(int vh) { + Point const& point = alpha_complex_.get_point(vh); + return pt_cgal_to_cython(point); + } + + bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, + bool default_filtration_value){ + return alpha_complex_.create_complex(*simplex_tree, max_alpha_square); + } + + private: + Alpha_complex_3d<Complexity, false, false> alpha_complex_; +}; + + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // INCLUDE_ALPHA_COMPLEX_FACTORY_H_ diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 3ac5db1f..38533a08 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -11,12 +11,8 @@ #ifndef INCLUDE_ALPHA_COMPLEX_INTERFACE_H_ #define INCLUDE_ALPHA_COMPLEX_INTERFACE_H_ -#include <gudhi/Simplex_tree.h> -#include <gudhi/Alpha_complex.h> -#include <CGAL/Epeck_d.h> -#include <CGAL/Epick_d.h> - -#include <boost/range/adaptor/transformed.hpp> +#include "Alpha_complex_factory.h" +#include <gudhi/Alpha_complex_options.h> #include "Simplex_tree_interface.h" @@ -30,67 +26,51 @@ namespace Gudhi { namespace alpha_complex { class Alpha_complex_interface { - 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; - - 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; - } - - template <typename CgalPointType> - static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) { - return CgalPointType(vec.size(), vec.begin(), vec.end()); - } - 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); + Alpha_complex_interface(const std::vector<std::vector<double>>& points, bool fast_version, bool exact_version) + : points_(points.begin(), points.end()), + fast_version_(fast_version), + exact_version_(exact_version) { } std::vector<double> get_point(int vh) { - 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); - } + return alpha_ptr_->get_point(vh); } - void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool exact_version, + void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, 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); + 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); + } + } } private: + std::unique_ptr<Abstract_alpha_complex> alpha_ptr_; + std::vector<std::vector<double>> points_; bool fast_version_; - std::unique_ptr<Alpha_complex<Exact_kernel>> ac_exact_ptr_; - std::unique_ptr<Alpha_complex<Inexact_kernel>> ac_inexact_ptr_; + bool exact_version_; }; } // namespace alpha_complex diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py index 943ad2c4..a4ee260b 100755 --- a/src/python/test/test_alpha_complex.py +++ b/src/python/test/test_alpha_complex.py @@ -8,7 +8,7 @@ - YYYY/MM Author: Description of the modification """ -from gudhi import AlphaComplex, SimplexTree +import gudhi as gd import math import numpy as np import pytest @@ -25,17 +25,16 @@ __license__ = "MIT" def _empty_alpha(precision): - alpha_complex = AlphaComplex(points=[[0, 0]], precision = precision) + alpha_complex = gd.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') + for precision in ['fast', 'safe', 'exact']: + _empty_alpha(precision) def _infinite_alpha(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - alpha_complex = AlphaComplex(points=point_list, precision = precision) + alpha_complex = gd.AlphaComplex(points=point_list, precision = precision) assert alpha_complex.__is_defined() == True simplex_tree = alpha_complex.create_simplex_tree() @@ -84,13 +83,12 @@ def _infinite_alpha(precision): assert False def test_infinite_alpha(): - _infinite_alpha('fast') - _infinite_alpha('safe') - _infinite_alpha('exact') + for precision in ['fast', 'safe', 'exact']: + _infinite_alpha(precision) def _filtered_alpha(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - filtered_alpha = AlphaComplex(points=point_list, precision = precision) + filtered_alpha = gd.AlphaComplex(points=point_list, precision = precision) simplex_tree = filtered_alpha.create_simplex_tree(max_alpha_square=0.25) @@ -128,9 +126,8 @@ def _filtered_alpha(precision): assert simplex_tree.get_cofaces([0], 1) == [([0, 1], 0.25), ([0, 2], 0.25)] def test_filtered_alpha(): - _filtered_alpha('fast') - _filtered_alpha('safe') - _filtered_alpha('exact') + for precision in ['fast', 'safe', 'exact']: + _filtered_alpha(precision) def _safe_alpha_persistence_comparison(precision): #generate periodic signal @@ -144,10 +141,10 @@ def _safe_alpha_persistence_comparison(precision): embedding2 = [[signal[i], delayed[i]] for i in range(len(time))] #build alpha complex and simplex tree - alpha_complex1 = AlphaComplex(points=embedding1, precision = precision) + alpha_complex1 = gd.AlphaComplex(points=embedding1, precision = precision) simplex_tree1 = alpha_complex1.create_simplex_tree() - alpha_complex2 = AlphaComplex(points=embedding2, precision = precision) + alpha_complex2 = gd.AlphaComplex(points=embedding2, precision = precision) simplex_tree2 = alpha_complex2.create_simplex_tree() diag1 = simplex_tree1.persistence() @@ -165,7 +162,7 @@ def test_safe_alpha_persistence_comparison(): def _delaunay_complex(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - filtered_alpha = AlphaComplex(points=point_list, precision = precision) + filtered_alpha = gd.AlphaComplex(points=point_list, precision = precision) simplex_tree = filtered_alpha.create_simplex_tree(default_filtration_value = True) @@ -197,6 +194,54 @@ def _delaunay_complex(precision): assert math.isnan(filtered_value[1]) def test_delaunay_complex(): - _delaunay_complex('fast') - _delaunay_complex('safe') - _delaunay_complex('exact') + for precision in ['fast', 'safe', 'exact']: + _delaunay_complex(precision) + +def _3d_points_on_a_plane(precision, default_filtration_value): + alpha = gd.AlphaComplex(off_file=gd.__root_source_dir__ + '/data/points/alphacomplexdoc.off', + precision = precision) + + simplex_tree = alpha.create_simplex_tree(default_filtration_value = default_filtration_value) + assert simplex_tree.dimension() == 2 + assert simplex_tree.num_vertices() == 7 + assert simplex_tree.num_simplices() == 25 + +def test_3d_points_on_a_plane(): + 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) + st_alpha = alpha.create_simplex_tree(default_filtration_value = False) + # New AlphaComplex for get_point to work + delaunay = gd.AlphaComplex(points=points, precision = precision) + st_delaunay = delaunay.create_simplex_tree(default_filtration_value = True) + + delaunay_tetra = [] + for sk in st_delaunay.get_skeleton(4): + if len(sk[0]) == 4: + tetra = [delaunay.get_point(sk[0][0]), + delaunay.get_point(sk[0][1]), + delaunay.get_point(sk[0][2]), + delaunay.get_point(sk[0][3]) ] + delaunay_tetra.append(sorted(tetra, key=lambda tup: tup[0])) + + alpha_tetra = [] + for sk in st_alpha.get_skeleton(4): + if len(sk[0]) == 4: + tetra = [alpha.get_point(sk[0][0]), + alpha.get_point(sk[0][1]), + alpha.get_point(sk[0][2]), + alpha.get_point(sk[0][3]) ] + alpha_tetra.append(sorted(tetra, key=lambda tup: tup[0])) + + # Check the tetrahedrons from one list are in the second one + assert len(alpha_tetra) == len(delaunay_tetra) + for tetra_from_del in delaunay_tetra: + assert tetra_from_del in alpha_tetra + +def test_3d_tetrahedrons(): + for precision in ['fast', 'safe', 'exact']: + _3d_tetrahedrons(precision) |