Factory to build python alpha complex. 3d version of alpha complex. Needs change in doc as alphacomplexdoc.off is a fake 3d off data file

4 files changed, 185 insertions, 107 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..7cab73bd
--- /dev/null
+++ b/src/python/include/Alpha_complex_factory.h
@@ -0,0 +1,144 @@
+/*    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& 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());
+}
+
+class Abstract_alpha_complex {
+ public:
+  virtual std::vector<double> get_point(int vh) = 0;
+  virtual void 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 Exact_kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+  using Point_exact_kernel = typename Exact_kernel::Point_d;
+
+ public:
+  Exact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version)
+      : exact_version_(exact_version) {
+    ac_exact_ptr_ = std::make_unique<Alpha_complex<Exact_kernel>>(
+      boost::adaptors::transform(points, pt_cython_to_cgal<Point_exact_kernel>));
+  }
+
+  std::vector<double> get_point(int vh) {
+    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,
+                           bool default_filtration_value){
+    ac_exact_ptr_->create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value);
+  }
+
+ private:
+  bool exact_version_;
+  std::unique_ptr<Alpha_complex<Exact_kernel>> ac_exact_ptr_;
+};
+
+class Inexact_Alphacomplex_dD : public Abstract_alpha_complex {
+ private:
+  using Inexact_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
+  using Point_inexact_kernel = typename Inexact_kernel::Point_d;
+
+ public:
+  Inexact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version)
+      : exact_version_(exact_version) {
+    ac_inexact_ptr_ = std::make_unique<Alpha_complex<Inexact_kernel>>(
+      boost::adaptors::transform(points, pt_cython_to_cgal<Point_inexact_kernel>));
+  }
+
+  std::vector<double> get_point(int vh) {
+    Point_inexact_kernel const& point = ac_inexact_ptr_->get_point(vh);
+    return pt_cgal_to_cython(point);
+  }
+  void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                           bool default_filtration_value){
+    ac_inexact_ptr_->create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value);
+  }
+
+ private:
+  bool exact_version_;
+  std::unique_ptr<Alpha_complex<Inexact_kernel>> ac_inexact_ptr_;
+};
+
+template <complexity Complexity>
+class Alphacomplex_3D : public Abstract_alpha_complex {
+ private:
+  using Point_3 = typename Alpha_complex_3d<Complexity, false, false>::Bare_point_3;
+
+  static Point_3 pt_cython_to_cgal_3(std::vector<double> const& vec) {
+    return Point_3(vec[0], vec[1], vec[2]);
+  }
+
+ public:
+  Alphacomplex_3D(const std::vector<std::vector<double>>& points) {
+    alpha3d_ptr_ = std::make_unique<Alpha_complex_3d<Complexity, false, false>>(
+      boost::adaptors::transform(points, pt_cython_to_cgal_3));
+  }
+
+  std::vector<double> get_point(int vh) {
+    Point_3 const& point = alpha3d_ptr_->get_point(vh);
+    return pt_cgal_to_cython(point);
+  }
+
+  void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                           bool default_filtration_value){
+    alpha3d_ptr_->create_complex(*simplex_tree, max_alpha_square);
+    if (default_filtration_value) {
+      // TODO
+    }
+  }
+
+ private:
+  std::unique_ptr<Alpha_complex_3d<Complexity, false, false>> alpha3d_ptr_;
+};
+
+
+}  // 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..c4b9df4b 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,35 @@ 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>));
+  Alpha_complex_interface(const std::vector<std::vector<double>>& points, bool fast_version, bool exact_version) {
+    if (points[0].size() == 3) {
+      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);
     } else {
-      ac_exact_ptr_ = std::make_unique<Alpha_complex<Exact_kernel>>(
-          boost::adaptors::transform(points, pt_cython_to_cgal<Point_exact_kernel>));
+      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_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) {
-    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);
+      alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value);
   }
 
  private:
-  bool fast_version_;
-  std::unique_ptr<Alpha_complex<Exact_kernel>> ac_exact_ptr_;
-  std::unique_ptr<Alpha_complex<Inexact_kernel>> ac_inexact_ptr_;
+  std::unique_ptr<Abstract_alpha_complex> alpha_ptr_;
 };
 
 }  // namespace alpha_complex