Add weight support

4 files changed, 194 insertions, 51 deletions

diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst
index b116944c..6f35cc15 100644
--- a/src/python/doc/alpha_complex_user.rst
+++ b/src/python/doc/alpha_complex_user.rst
@@ -44,16 +44,15 @@ This example builds the alpha-complex from the given points:
 
 .. testcode::
 
-    import gudhi
-    alpha_complex = gudhi.AlphaComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]])
-
-    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)
+    from gudhi import AlphaComplex
+    ac = AlphaComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]])
+
+    stree = ac.create_simplex_tree()
+    print('Alpha complex is of dimension ', stree.dimension(), ' - ',
+      stree.num_simplices(), ' simplices - ', stree.num_vertices(), ' vertices.')
+
     fmt = '%s -> %.2f'
-    for filtered_value in simplex_tree.get_filtration():
+    for filtered_value in stree.get_filtration():
         print(fmt % tuple(filtered_value))
 
 The output is:
@@ -174,6 +173,78 @@ of speed-up, since we always first build the full filtered complex, so it is rec
 :paramref:`~gudhi.AlphaComplex.create_simplex_tree.max_alpha_square`.
 In the following example, a threshold of :math:`\alpha^2 = 32.0` is used.
 
+Weighted specific version
+^^^^^^^^^^^^^^^^^^^^^^^^^
+
+:Requires: `Eigen <installation.html#eigen>`_ :math:`\geq` 3.1.0 and `CGAL <installation.html#cgal>`_ :math:`\geq` 5.1.0.
+
+A weighted version for Alpha complex is available. It is like a usual Alpha complex, but based on a
+`CGAL regular triangulation <https://doc.cgal.org/latest/Triangulation/index.html#title20>`_
+of Delaunay.
+
+This example builds the CGAL weighted alpha shapes from a small molecule, and initializes the alpha complex with
+it. This example is taken from
+`CGAL 3d weighted alpha shapes <https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title13>`_.
+
+Then, it is asked to display information about the alpha complex.
+
+.. testcode::
+
+    from gudhi import AlphaComplex
+    wgt_ac = AlphaComplex(weighted_points=[[[ 1., -1., -1.], 4.],
+                                           [[-1.,  1., -1.], 4.],
+                                           [[-1., -1.,  1.], 4.],
+                                           [[ 1.,  1.,  1.], 4.],
+                                           [[ 2.,  2.,  2.], 1.]])
+    # equivalent to:
+    # wgt_ac = AlphaComplex(points=[[ 1., -1., -1.],
+    #                               [-1.,  1., -1.],
+    #                               [-1., -1.,  1.],
+    #                               [ 1.,  1.,  1.],
+    #                               [ 2.,  2.,  2.]],
+    #                       weights = [4., 4., 4., 4., 1.])
+
+    stree = wgt_ac.create_simplex_tree()
+    print('Weighted alpha complex is of dimension ', stree.dimension(), ' - ',
+          stree.num_simplices(), ' simplices - ', stree.num_vertices(), ' vertices.')
+    fmt = '%s -> %.2f'
+    for filtered_value in stree.get_filtration():
+        print(fmt % tuple(filtered_value))
+
+The output is:
+
+.. testoutput::
+
+   Weighted alpha complex is of dimension 3 - 29 simplices - 5 vertices.
+   [ 0 ] -> [-4] 
+   [ 1 ] -> [-4] 
+   [ 2 ] -> [-4] 
+   [ 3 ] -> [-4] 
+   [ 1, 0 ] -> [-2] 
+   [ 2, 0 ] -> [-2] 
+   [ 2, 1 ] -> [-2] 
+   [ 3, 0 ] -> [-2] 
+   [ 3, 1 ] -> [-2] 
+   [ 3, 2 ] -> [-2] 
+   [ 2, 1, 0 ] -> [-1.33333] 
+   [ 3, 1, 0 ] -> [-1.33333] 
+   [ 3, 2, 0 ] -> [-1.33333] 
+   [ 3, 2, 1 ] -> [-1.33333] 
+   [ 3, 2, 1, 0 ] -> [-1] 
+   [ 4 ] -> [-1] 
+   [ 4, 2 ] -> [-1] 
+   [ 4, 0 ] -> [23] 
+   [ 4, 1 ] -> [23] 
+   [ 4, 2, 0 ] -> [23] 
+   [ 4, 2, 1 ] -> [23] 
+   [ 4, 3 ] -> [23] 
+   [ 4, 3, 2 ] -> [23] 
+   [ 4, 1, 0 ] -> [95] 
+   [ 4, 2, 1, 0 ] -> [95] 
+   [ 4, 3, 0 ] -> [95] 
+   [ 4, 3, 1 ] -> [95] 
+   [ 4, 3, 2, 0 ] -> [95] 
+   [ 4, 3, 2, 1 ] -> [95] 
 
 Example from OFF file
 ^^^^^^^^^^^^^^^^^^^^^
@@ -186,14 +257,9 @@ Then, it computes the persistence diagram and displays it:
    :include-source:
 
     import matplotlib.pyplot as plt
-    import gudhi
-    alpha_complex = gudhi.AlphaComplex(off_file=gudhi.__root_source_dir__ + \
-        '/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)
-    diag = simplex_tree.persistence()
-    gudhi.plot_persistence_diagram(diag)
+    import gudhi as gd
+    ac = gd.AlphaComplex(off_file=gd.__root_source_dir__ + '/data/points/tore3D_300.off')
+    stree = ac.create_simplex_tree()
+    dgm = stree.persistence()
+    gd.plot_persistence_diagram(dgm, legend = True)
     plt.show()
diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx
index b7c20f74..681faebe 100644
--- a/src/python/gudhi/alpha_complex.pyx
+++ b/src/python/gudhi/alpha_complex.pyx
@@ -29,7 +29,7 @@ __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, bool exact_version) nogil except +
+        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 +
 
@@ -56,26 +56,35 @@ cdef class AlphaComplex:
     cdef Alpha_complex_interface * this_ptr
 
     # Fake constructor that does nothing but documenting the constructor
-    def __init__(self, points=[], off_file='', weights=[], weight_file='', precision='safe'):
+    def __init__(self, points=[], off_file='', weights=[], weight_file='', weighted_points=[],
+                 precision='safe'):
         """AlphaComplex constructor.
 
         :param points: A list of points in d-Dimension.
-        :type points: list of list of double
+        :type points: Iterable[Iterable[float]]
 
-        :param off_file: An `OFF file style <fileformats.html#off-file-format>`_ name. `points` are
-            read and overwritten by the points in the `off_file`.
+        :param off_file: An `OFF file style <fileformats.html#off-file-format>`_ name. 
+            If an `off_file` is given with `points` or `weighted_points`, only points from the
+            file are taken into account.
         :type off_file: string
 
         :param weights: A list of weights. If set, the number of weights must correspond to the
             number of points.
-        :type weights: list of double
+        :type weights: Iterable[float]
 
         :param weight_file: A file containing a list of weights (one per line).
-            `weights` are read and overwritten by the weights in the `weight_file`.
-            If set, the number of weights must correspond to the number of points.
+            If a `weight_file` is given with `weights` or `weighted_points`, only weights from the
+            file are taken into account.
+
         :type weight_file: string
 
-        :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'.
+        :param weighted_points: A list of points in d-Dimension and its weight.
+            If `weighted_points` are given with `weights` or `points`, these last ones will
+            not be taken into account.
+        :type weighted_points: Iterable[Iterable[float], float]
+
+        :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is
+            'safe'.
         :type precision: string
 
         :raises FileNotFoundError: If `off_file` and/or `weight_file` is set but not found.
@@ -83,11 +92,16 @@ cdef class AlphaComplex:
         """
 
     # The real cython constructor
-    def __cinit__(self, points = [], off_file = '', weights=[], weight_file='', precision = 'safe'):
+    def __cinit__(self, points = [], off_file = '', weights=[], weight_file='', weighted_points=[],
+                  precision = 'safe'):
         assert precision in ['fast', 'safe', 'exact'], "Alpha complex precision can only be 'fast', 'safe' or 'exact'"
         cdef bool fast = precision == 'fast'
         cdef bool exact = precision == 'exact'
 
+        if len(weighted_points) > 0:
+            points = [wpt[0] for wpt in weighted_points]
+            weights = [wpt[1] for wpt in weighted_points]
+
         if off_file:
             if os.path.isfile(off_file):
                 points = read_points_from_off_file(off_file = off_file)
@@ -100,20 +114,18 @@ cdef class AlphaComplex:
             else:
                 raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), weight_file)
 
+        # weights are set but is inconsistent with the number of points
+        if len(weights) != 0 and len(weights) != len(points):
+            raise ValueError("Inconsistency between the number of points and weights")
+
         # need to copy the points to use them without the gil
         cdef vector[vector[double]] pts
         cdef vector[double] wgts
         pts = points
+        wgts = weights
         if len(weights) == 0:
             with nogil:
-                self.this_ptr = new Alpha_complex_interface(pts, fast, exact)
-        else:
-            if len(weights) == len(points):
-                wgts = weights
-                with nogil:
-                    self.this_ptr = new Alpha_complex_interface(pts, fast, exact)
-            else:
-                raise ValueError("Inconsistency between the number of points and weights")
+                self.this_ptr = new Alpha_complex_interface(pts, wgts, fast, exact)
 
     def __dealloc__(self):
         if self.this_ptr != NULL:
diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h
index 3405fdd6..36e98615 100644
--- a/src/python/include/Alpha_complex_factory.h
+++ b/src/python/include/Alpha_complex_factory.h
@@ -55,13 +55,13 @@ class Abstract_alpha_complex {
   virtual ~Abstract_alpha_complex() = default;
 };
 
-class Exact_Alphacomplex_dD final : public Abstract_alpha_complex {
+class Exact_alpha_complex_dD final : public Abstract_alpha_complex {
  private:
   using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
   using Point = typename Kernel::Point_d;
 
  public:
-  Exact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version)
+  Exact_alpha_complex_dD(const std::vector<std::vector<double>>& points, bool exact_version)
     : exact_version_(exact_version),
       alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>)) {
   }
@@ -81,13 +81,13 @@ class Exact_Alphacomplex_dD final : public Abstract_alpha_complex {
   Alpha_complex<Kernel> alpha_complex_;
 };
 
-class Inexact_Alphacomplex_dD final : public Abstract_alpha_complex {
+class Inexact_alpha_complex_dD final : public Abstract_alpha_complex {
  private:
   using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
   using Point = typename Kernel::Point_d;
 
  public:
-  Inexact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version)
+  Inexact_alpha_complex_dD(const std::vector<std::vector<double>>& points, bool exact_version)
     : exact_version_(exact_version),
       alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>)) {
   }
@@ -106,8 +106,61 @@ class Inexact_Alphacomplex_dD final : public Abstract_alpha_complex {
   Alpha_complex<Kernel> alpha_complex_;
 };
 
+class Exact_weighted_alpha_complex_dD final : public Abstract_alpha_complex {
+ private:
+  using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+  using Point = typename Kernel::Point_d;
+
+ public:
+  Exact_weighted_alpha_complex_dD(const std::vector<std::vector<double>>& points,
+                                  const std::vector<double>& weights, bool exact_version)
+    : exact_version_(exact_version),
+      alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>), weights) {
+  }
+
+  virtual std::vector<double> get_point(int vh) override {
+    Point const& point = alpha_complex_.get_point(vh).point();
+    return pt_cgal_to_cython(point);
+  }
+
+  virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                                   bool default_filtration_value) override {
+    return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value);
+  }
+
+ private:
+  bool exact_version_;
+  Alpha_complex<Kernel, true> alpha_complex_;
+};
+
+class Inexact_weighted_alpha_complex_dD final : public Abstract_alpha_complex {
+ private:
+  using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
+  using Point = typename Kernel::Point_d;
+
+ public:
+  Inexact_weighted_alpha_complex_dD(const std::vector<std::vector<double>>& points,
+                                    const std::vector<double>& weights, bool exact_version)
+    : exact_version_(exact_version),
+      alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>), weights) {
+  }
+
+  virtual std::vector<double> get_point(int vh) override {
+    Point const& point = alpha_complex_.get_point(vh).point();
+    return pt_cgal_to_cython(point);
+  }
+  virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square,
+                                   bool default_filtration_value) override {
+    return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value);
+  }
+
+ private:
+  bool exact_version_;
+  Alpha_complex<Kernel, true> alpha_complex_;
+};
+
 template <complexity Complexity>
-class Alphacomplex_3D final : public Abstract_alpha_complex {
+class Alpha_complex_3D final : public Abstract_alpha_complex {
  private:
   using Point = typename Alpha_complex_3d<Complexity, false, false>::Bare_point_3;
 
@@ -116,7 +169,7 @@ class Alphacomplex_3D final : public Abstract_alpha_complex {
   }
 
  public:
-  Alphacomplex_3D(const std::vector<std::vector<double>>& points)
+  Alpha_complex_3D(const std::vector<std::vector<double>>& points)
     : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3)) {
   }
 
diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h
index 23be194d..43c96b2f 100644
--- a/src/python/include/Alpha_complex_interface.h
+++ b/src/python/include/Alpha_complex_interface.h
@@ -27,8 +27,11 @@ namespace alpha_complex {
 
 class Alpha_complex_interface {
  public:
-  Alpha_complex_interface(const std::vector<std::vector<double>>& points, bool fast_version, bool exact_version)
+  Alpha_complex_interface(const std::vector<std::vector<double>>& points,
+                          const std::vector<double>& weights,
+                          bool fast_version, bool exact_version)
   : points_(points),
+    weights_(weights),
     fast_version_(fast_version),
     exact_version_(exact_version) {
   }
@@ -41,13 +44,13 @@ class Alpha_complex_interface {
                            bool default_filtration_value) {
     if (points_.size() > 0) {
       std::size_t dimension = points_[0].size();
-      if (dimension == 3 && !default_filtration_value) {
+      if (dimension == 3 && weights_.size() == 0 && !default_filtration_value) {
         if (fast_version_)
-          alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::FAST>>(points_);
+          alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::FAST>>(points_);
         else if (exact_version_)
-          alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::EXACT>>(points_);
+          alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::EXACT>>(points_);
         else
-          alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::SAFE>>(points_);
+          alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::SAFE>>(points_);
         if (!alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value)) {
           // create_simplex_tree will fail if all points are on a plane - Retry with dD by setting dimension to 2
           dimension--;
@@ -55,11 +58,19 @@ class Alpha_complex_interface {
         }
       }
       // Not ** else ** because we have to take into account if 3d fails
-      if (dimension != 3 || default_filtration_value) {
+      if (dimension != 3 || weights_.size() != 0 || default_filtration_value) {
         if (fast_version_) {
-          alpha_ptr_ = std::make_unique<Inexact_Alphacomplex_dD>(points_, exact_version_);
+          if (weights_.size() == 0) {
+            alpha_ptr_ = std::make_unique<Inexact_alpha_complex_dD>(points_, exact_version_);
+          } else {
+            alpha_ptr_ = std::make_unique<Inexact_weighted_alpha_complex_dD>(points_, weights_, exact_version_);
+          }
         } else {
-          alpha_ptr_ = std::make_unique<Exact_Alphacomplex_dD>(points_, exact_version_);
+          if (weights_.size() == 0) {
+            alpha_ptr_ = std::make_unique<Exact_alpha_complex_dD>(points_, exact_version_);
+          } else {
+            alpha_ptr_ = std::make_unique<Exact_weighted_alpha_complex_dD>(points_, weights_, exact_version_);
+          }
         }
         alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value);
       }
@@ -69,6 +80,7 @@ class Alpha_complex_interface {
  private:
   std::unique_ptr<Abstract_alpha_complex> alpha_ptr_;
   std::vector<std::vector<double>> points_;
+  std::vector<double> weights_;
   bool fast_version_;
   bool exact_version_;
 };