Merge pull request #169 from VincentRouvreau/exact_alpha_complex_python_version

Fix #107 - Exact alpha complex python version

4 files changed, 48 insertions, 11 deletions

diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 9af85eac..b558d4c4 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -114,9 +114,9 @@ if(PYTHONINTERP_FOUND)
       set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'nerve_gic', ")
     endif ()
     if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
+      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'alpha_complex', ")
       set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'subsampling', ")
       set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'tangential_complex', ")
-      set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'alpha_complex', ")
       set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'euclidean_witness_complex', ")
       set(GUDHI_PYTHON_MODULES_TO_COMPILE "${GUDHI_PYTHON_MODULES_TO_COMPILE}'euclidean_strong_witness_complex', ")
     endif ()
@@ -162,10 +162,21 @@ if(PYTHONINTERP_FOUND)
           set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_GMPXX', ")
           add_GUDHI_PYTHON_lib("${GMPXX_LIBRARIES}")
           set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${GMPXX_LIBRARIES_DIR}', ")
-        message("** Add gmpxx ${GMPXX_LIBRARIES_DIR}")
+          message("** Add gmpxx ${GMPXX_LIBRARIES_DIR}")
         endif(GMPXX_FOUND)
       endif(GMP_FOUND)
-    endif(CGAL_FOUND)
+      if(MPFR_FOUND)
+        add_gudhi_debug_info("MPFR_LIBRARIES = ${MPFR_LIBRARIES}")
+        set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_MPFR', ")
+        add_GUDHI_PYTHON_lib("${MPFR_LIBRARIES}")
+        # In case CGAL is not header only, all MPFR variables are set except MPFR_LIBRARIES_DIR - Just set it
+        if(NOT MPFR_LIBRARIES_DIR)
+          get_filename_component(MPFR_LIBRARIES_DIR ${MPFR_LIBRARIES} PATH)
+        endif(NOT MPFR_LIBRARIES_DIR)
+        set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${MPFR_LIBRARIES_DIR}', ")
+        message("** Add mpfr ${MPFR_LIBRARIES}")
+      endif(MPFR_FOUND)
+endif(CGAL_FOUND)
 
     # Specific for Mac
     if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
@@ -229,7 +240,6 @@ if(PYTHONINTERP_FOUND)
                COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}"
                ${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)
-
       if(MATPLOTLIB_FOUND AND NUMPY_FOUND)
         # Tangential
         add_test(NAME tangential_complex_plain_homology_from_off_file_example_py_test
@@ -305,7 +315,6 @@ if(PYTHONINTERP_FOUND)
                WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
                COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}"
                ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/alpha_complex_from_points_example.py")
-
       if(MATPLOTLIB_FOUND AND NUMPY_FOUND)
         add_test(NAME alpha_complex_diagram_persistence_from_off_file_example_py_test
                WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
@@ -313,9 +322,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 -a 0.6)
       endif()
-
       add_gudhi_py_test(test_alpha_complex)
-
     endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
 
     if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx
index 8f2c98d5..24e36bea 100644
--- a/src/python/gudhi/alpha_complex.pyx
+++ b/src/python/gudhi/alpha_complex.pyx
@@ -66,7 +66,7 @@ cdef class AlphaComplex:
         """
 
     # The real cython constructor
-    def __cinit__(self, points=None, off_file=''):
+    def __cinit__(self, points = None, off_file = ''):
         if off_file:
             if os.path.isfile(off_file):
                 self.thisptr = new Alpha_complex_interface(str.encode(off_file), True)
@@ -99,7 +99,7 @@ cdef class AlphaComplex:
         cdef vector[double] point = self.thisptr.get_point(vertex)
         return point
 
-    def create_simplex_tree(self, max_alpha_square=float('inf')):
+    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
diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h
index 96353cc4..e9bbadb0 100644
--- a/src/python/include/Alpha_complex_interface.h
+++ b/src/python/include/Alpha_complex_interface.h
@@ -13,6 +13,7 @@
 
 #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>
@@ -28,7 +29,7 @@ namespace Gudhi {
 namespace alpha_complex {
 
 class Alpha_complex_interface {
-  using Dynamic_kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >;
+  using Dynamic_kernel = CGAL::Epeck_d< CGAL::Dynamic_dimension_tag >;
   using Point_d = Dynamic_kernel::Point_d;
 
  public:
@@ -51,7 +52,7 @@ class Alpha_complex_interface {
     std::vector<double> vd;
     try {
       Point_d const& ph = alpha_complex_->get_point(vh);
-      for (auto coord = ph.cartesian_begin(); coord < ph.cartesian_end(); coord++)
+      for (auto coord = ph.cartesian_begin(); coord != ph.cartesian_end(); coord++)
         vd.push_back(CGAL::to_double(*coord));
     } catch (std::out_of_range const&) {
       // std::out_of_range is thrown in case not found. Other exceptions must be re-thrown
diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py
index 24f8bf53..9b27fff2 100755
--- a/src/python/test/test_alpha_complex.py
+++ b/src/python/test/test_alpha_complex.py
@@ -1,4 +1,8 @@
 from gudhi import AlphaComplex, SimplexTree
+import math
+import numpy as np
+import itertools
+import pytest
 
 """ 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.
@@ -88,3 +92,28 @@ 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():
+    #generate periodic signal
+    time = np.arange(0, 10, 1)
+    signal = [math.sin(x) for x in time]
+    delta = math.pi
+    delayed = [math.sin(x + delta) for x in time]
+    
+    #construct embedding
+    embedding1 = [[signal[i], -signal[i]] for i in range(len(time))]
+    embedding2 = [[signal[i], delayed[i]] for i in range(len(time))]
+    
+    #build alpha complex and simplex tree
+    alpha_complex1 = AlphaComplex(points=embedding1)
+    simplex_tree1 = alpha_complex1.create_simplex_tree()
+    
+    alpha_complex2 = AlphaComplex(points=embedding2)
+    simplex_tree2 = alpha_complex2.create_simplex_tree()
+    
+    diag1 = simplex_tree1.persistence()
+    diag2 = simplex_tree2.persistence()
+
+    for (first_p, second_p) in itertools.zip_longest(diag1, diag2):
+        assert first_p[0] == pytest.approx(second_p[0])
+        assert first_p[1] == pytest.approx(second_p[1])