Merge branch 'master' into stree_iterator_boundary_and_opposite_vertex

7 files changed, 88 insertions, 4 deletions

diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
index df5c630e..b3dbc9bb 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
@@ -57,8 +57,6 @@ namespace Gudhi {
 
 namespace alpha_complex {
 
-thread_local double RELATIVE_PRECISION_OF_TO_DOUBLE = 0.00001;
-
 // Value_from_iterator returns the filtration value from an iterator on alpha shapes values
 //
 // FAST                         SAFE                         EXACT
diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake
index 7c982b3b..6a94d1f5 100644
--- a/src/cmake/modules/GUDHI_third_party_libraries.cmake
+++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake
@@ -19,6 +19,15 @@ if(GMP_FOUND)
   endif()
 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)
+  message("++ User explicit demand to force EIGEN_DEFAULT_DENSE_INDEX_TYPE to int")
+  add_definitions(-DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int)
+endif()
+
 # In CMakeLists.txt, when include(${CGAL_USE_FILE}), CMAKE_CXX_FLAGS are overwritten.
 # cf. http://doc.cgal.org/latest/Manual/installation.html#title40
 # A workaround is to include(${CGAL_USE_FILE}) before adding "-std=c++11".
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/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/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)