Merge pull request #326 from VincentRouvreau/alpha_complex_python_fast_and_exact

Alpha complex python enhancement

5 files changed, 163 insertions, 51 deletions

diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
index 7c898dfd..e17831d9 100644
--- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
@@ -11,6 +11,7 @@
 #include <boost/program_options.hpp>
 
 #include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h>
 
 #include <gudhi/Alpha_complex.h>
 #include <gudhi/Persistent_cohomology.h>
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')