diff options
author | Vincent Rouvreau <vincent.rouvreau@inria.fr> | 2022-04-22 14:40:21 +0200 |
---|---|---|
committer | Vincent Rouvreau <vincent.rouvreau@inria.fr> | 2022-04-22 14:40:21 +0200 |
commit | f2e139123b15fee6e1512aac3824dd6664561cdb (patch) | |
tree | 1f0a7b458a461a75c9f7fe3c43ee6980adf10335 /src/python | |
parent | 17dc48527dcc8ee7e5eab95f9fdde3e236f4ad47 (diff) | |
parent | 5857c571388a4349934a266ca187b2c2ac10818c (diff) |
Merge branch 'master' into simplex_tree_extended_persistence_enhancement
Diffstat (limited to 'src/python')
-rw-r--r-- | src/python/CMakeLists.txt | 12 | ||||
-rw-r--r-- | src/python/doc/alpha_complex_user.rst | 3 | ||||
-rw-r--r-- | src/python/gudhi/alpha_complex.pyx | 29 | ||||
-rw-r--r-- | src/python/gudhi/simplex_tree.pxd | 3 | ||||
-rw-r--r-- | src/python/gudhi/simplex_tree.pyx | 24 | ||||
-rw-r--r-- | src/python/include/Alpha_complex_interface.h | 10 | ||||
-rw-r--r-- | src/python/include/Simplex_tree_interface.h | 8 | ||||
-rwxr-xr-x | src/python/test/test_alpha_complex.py | 27 | ||||
-rwxr-xr-x | src/python/test/test_simplex_tree.py | 45 |
9 files changed, 159 insertions, 2 deletions
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index fbfb0d1f..e31af02e 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -180,6 +180,15 @@ 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) @@ -215,13 +224,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}") diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst index cfd22742..b060c86e 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 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/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index a8ed6d50..5642f82d 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -75,6 +75,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_interface<Gudhi::Simplex_tree_options_full_featured>>": diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 1bbf1539..97c7d2c8 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 @@ -473,6 +476,27 @@ cdef class SimplexTree: self.pcohptr.compute_persistence(homology_coeff_field, -1.) return self.pcohptr.compute_extended_persistence_subdiagrams(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. 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 dc9d01d7..aec29f9b 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -42,6 +42,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: @@ -165,6 +166,13 @@ class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> { #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 Complex_simplex_iterator get_simplices_iterator_begin() { // this specific case works because the range is just a pair of iterators - won't work if range was a vector 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_simplex_tree.py b/src/python/test/test_simplex_tree.py index 23458eb2..7abb4366 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -519,3 +519,48 @@ def test_simplex_tree_deep_copy_constructor(): 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. |