Merge pull request #393 from VincentRouvreau/weighted_alpha_complex_dD

Weighted alpha complex dD

17 files changed, 815 insertions, 117 deletions

diff --git a/.github/next_release.md b/.github/next_release.md
index cd2376eb..190f8408 100644
--- a/.github/next_release.md
+++ b/.github/next_release.md
@@ -1,13 +1,13 @@
 We are pleased to announce the release 3.4.0 of the GUDHI library.
 
-As a major new feature, the GUDHI library now offers ...
+As a major new feature, the GUDHI library now offers dD weighted alpha complex, ...
 
 We are now using GitHub to develop the GUDHI library, do not hesitate to [fork the GUDHI project on GitHub](https://github.com/GUDHI/gudhi-devel). From a user point of view, we recommend to download GUDHI user version (gudhi.3.4.0.tar.gz).
 
 Below is a list of changes made since GUDHI 3.3.0:
 
-- [Module](link)
-     - ...
+- [Alpha complex](https://gudhi.inria.fr/doc/latest/group__alpha__complex.html)
+     - the C++ weighted version for alpha complex is now available in dimension D.
 
 - [Module](link)
      - ...
diff --git a/Dockerfile_for_circleci_image b/Dockerfile_for_circleci_image
index 87f57071..ec1b8ff8 100644
--- a/Dockerfile_for_circleci_image
+++ b/Dockerfile_for_circleci_image
@@ -41,7 +41,6 @@ RUN apt-get install -y make \
     libgmp3-dev \
     libmpfr-dev \
     libtbb-dev \
-    libcgal-dev \
     locales \
     python3 \
     python3-pip \
@@ -51,6 +50,15 @@ RUN apt-get install -y make \
     pkg-config \
     curl
 
+RUN curl -LO "https://github.com/CGAL/cgal/releases/download/v5.1/CGAL-5.1.tar.xz" \
+    && tar xf CGAL-5.1.tar.xz \
+    && mkdir build \
+    && cd build \
+    && cmake -DCMAKE_BUILD_TYPE=Release ../CGAL-5.1/ \
+    && make install \
+    && cd .. \
+    && rm -rf build CGAL-5.1
+
 ADD .github/build-requirements.txt /
 ADD .github/test-requirements.txt /
 
diff --git a/Dockerfile_gudhi_installation b/Dockerfile_gudhi_installation
index 92430fce..ebd21f8d 100644
--- a/Dockerfile_gudhi_installation
+++ b/Dockerfile_gudhi_installation
@@ -23,6 +23,9 @@ ENV LANG en_US.UTF-8
 ENV LANGUAGE en_US:en
 ENV LC_ALL en_US.UTF-8
 
+# Update again
+RUN apt-get update
+
 # Required for Gudhi compilation
 RUN apt-get install -y make \
     g++ \
@@ -47,6 +50,15 @@ RUN apt-get install -y make \
     pkg-config \
     curl
 
+RUN curl -LO "https://github.com/CGAL/cgal/releases/download/v5.1/CGAL-5.1.tar.xz" \
+    && tar xf CGAL-5.1.tar.xz \
+    && mkdir build \
+    && cd build \
+    && cmake -DCMAKE_BUILD_TYPE=Release ../CGAL-5.1/ \
+    && make install \
+    && cd .. \
+    && rm -rf build CGAL-5.1
+
 RUN pip3 install \
     numpy \
     matplotlib \
diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h
index 60da7169..c068b268 100644
--- a/src/Alpha_complex/doc/Intro_alpha_complex.h
+++ b/src/Alpha_complex/doc/Intro_alpha_complex.h
@@ -22,6 +22,18 @@ namespace alpha_complex {
  *
  * @{
  *
+<div class="toc">
+Table of Contents
+<ul>
+<li class="level1"><a href="#definition">Definition</a></li>
+<li class="level1"><a href="#pointsexample">Example from points</a></li>
+<li class="level1"><a href="#createcomplexalgorithm">Create complex algorithm</a></li>
+<li class="level1"><a href="#weightedversion">Weighted specific version</a></li>
+<li class="level1"><a href="#offexample">Example from OFF file</a></li>
+<li class="level1"><a href="#weighted3dexample">3d specific version</a></li>
+</ul>
+</div>
+
  * \section definition Definition
  *
  * Alpha_complex is a <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a>
@@ -62,7 +74,7 @@ namespace alpha_complex {
  * [10^12,10^12+10^6]. Using `CGAL::Epick_d` 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 \ref cgal &ge; 5.0.0.
+ * - For performances reasons, it is advised to use \ref eigen &ge; 3.3.5 and \ref cgal &ge; 5.2.0.
  *
  * \section pointsexample Example from points
  *
@@ -146,6 +158,30 @@ namespace alpha_complex {
  * `SimplicialComplexForAlpha::prune_above_filtration()`).
  * In the following example, the value is given by the user as argument of the program.
  *
+ * \section weightedversion Weighted specific version
+ * <b>Requires:</b> \ref eigen &ge; 3.1.0 and \ref cgal &ge; 5.1.0.
+ * 
+ * A weighted version for Alpha complex is available (cf. Alpha_complex). It is like a usual Alpha complex, but based
+ * on a <a href="https://doc.cgal.org/latest/Triangulation/index.html#title20">CGAL regular triangulation</a> instead
+ * 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 <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title13">CGAL 3d
+ * weighted alpha shapes</a>.
+ *
+ * Then, it is asked to display information about the alpha complex.
+ *
+ * \include Alpha_complex/Weighted_alpha_complex_from_points.cpp
+ *
+ * When launching:
+ *
+ * \code $> ./Weighted_alpha_complex_example_from_points
+ * \endcode
+ *
+ * the program output is:
+ *
+ * \include Alpha_complex/weightedalpha3dfrompoints_for_doc.txt
+ *
  *
  * \section offexample Example from OFF file
  *
@@ -166,7 +202,7 @@ namespace alpha_complex {
  * \include Alpha_complex/alphaoffreader_for_doc_32.txt
  *
  *
- * \section weighted3dexample 3d specific example
+ * \section weighted3dexample 3d specific version
  *
  * A specific module for Alpha complex is available in 3d (cf. Alpha_complex_3d) and allows to construct standard,
  * weighted, periodic or weighted and periodic versions of alpha complexes. Alpha values computation can be
@@ -181,14 +217,7 @@ namespace alpha_complex {
  *
  * \include Alpha_complex/Weighted_alpha_complex_3d_from_points.cpp
  *
- * When launching:
- *
- * \code $> ./Alpha_complex_example_weighted_3d_from_points
- * \endcode
- *
- * the program output is:
- *
- * \include Alpha_complex/weightedalpha3dfrompoints_for_doc.txt
+ * The results will be the same as in \ref weightedversion .
  *
  */
 /** @} */  // end defgroup alpha_complex
diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt
index 6e231773..1fc2330a 100644
--- a/src/Alpha_complex/example/CMakeLists.txt
+++ b/src/Alpha_complex/example/CMakeLists.txt
@@ -44,7 +44,7 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
         ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt)
     set_tests_properties(Alpha_complex_example_fast_from_off_32_diff_files PROPERTIES DEPENDS Alpha_complex_example_fast_from_off_32) 
   endif()
- endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
+endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
 
 if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
   add_executable ( Alpha_complex_example_weighted_3d_from_points Weighted_alpha_complex_3d_from_points.cpp )
@@ -64,3 +64,12 @@ if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
           COMMAND $<TARGET_FILE:Alpha_complex_example_3d_from_points>)
 
 endif(NOT CGAL_VERSION VERSION_LESS 4.11.0)
+
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
+  add_executable ( Weighted_alpha_complex_example_from_points Weighted_alpha_complex_from_points.cpp )
+  target_link_libraries(Weighted_alpha_complex_example_from_points  ${CGAL_LIBRARY})
+  if (TBB_FOUND)
+    target_link_libraries(Weighted_alpha_complex_example_from_points ${TBB_LIBRARIES})
+  endif()
+  add_test(NAME Weighted_alpha_complex_example_from_points COMMAND $<TARGET_FILE:Weighted_alpha_complex_example_from_points>)
+endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
diff --git a/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp
index c044194e..ee12d418 100644
--- a/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp
+++ b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp
@@ -7,7 +7,7 @@
 #include <vector>
 #include <limits>  // for numeric limits
 
-// Complexity = FAST, weighted = true, periodic = false
+// Complexity = SAFE, weighted = true, periodic = false
 using Weighted_alpha_complex_3d =
     Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>;
 using Bare_point = Weighted_alpha_complex_3d::Bare_point_3;
@@ -18,11 +18,11 @@ int main(int argc, char **argv) {
   // Init of a list of points and weights from a small molecule
   // ----------------------------------------------------------------------------
   std::vector<Weighted_point> weighted_points;
-  weighted_points.push_back(Weighted_point(Bare_point(1, -1, -1), 4.));
-  weighted_points.push_back(Weighted_point(Bare_point(-1, 1, -1), 4.));
-  weighted_points.push_back(Weighted_point(Bare_point(-1, -1, 1), 4.));
-  weighted_points.push_back(Weighted_point(Bare_point(1, 1, 1), 4.));
-  weighted_points.push_back(Weighted_point(Bare_point(2, 2, 2), 1.));
+  weighted_points.emplace_back(Bare_point(1, -1, -1), 4.);
+  weighted_points.emplace_back(Bare_point(-1, 1, -1), 4.);
+  weighted_points.emplace_back(Bare_point(-1, -1, 1), 4.);
+  weighted_points.emplace_back(Bare_point(1, 1, 1),   4.);
+  weighted_points.emplace_back(Bare_point(2, 2, 2),   1.);
 
   // ----------------------------------------------------------------------------
   // Init of an alpha complex from the list of points
@@ -34,10 +34,10 @@ int main(int argc, char **argv) {
     // ----------------------------------------------------------------------------
     // Display information about the alpha complex
     // ----------------------------------------------------------------------------
-    std::clog << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
+    std::clog << "Weighted alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
               << " simplices - " << simplex.num_vertices() << " vertices." << std::endl;
 
-    std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
+    std::clog << "Iterator on weighted alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
     for (auto f_simplex : simplex.filtration_simplex_range()) {
       std::clog << "   ( ";
       for (auto vertex : simplex.simplex_vertex_range(f_simplex)) {
diff --git a/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp b/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp
new file mode 100644
index 00000000..d1f3e436
--- /dev/null
+++ b/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp
@@ -0,0 +1,52 @@
+#include <gudhi/Alpha_complex.h>
+// to construct a simplex_tree from alpha complex
+#include <gudhi/Simplex_tree.h>
+
+#include <CGAL/Epeck_d.h>
+
+#include <iostream>
+#include <vector>
+
+// Explicit dimension 2 Epeck_d kernel
+using Kernel = CGAL::Epeck_d< CGAL::Dimension_tag<3> >;
+using Bare_point = Kernel::Point_d;
+using Weighted_point = Kernel::Weighted_point_d;
+using Vector_of_points = std::vector<Weighted_point>;
+
+int main() {
+  // ----------------------------------------------------------------------------
+  // Init of a list of points and weights from a small molecule
+  // ----------------------------------------------------------------------------
+  Vector_of_points points;
+  points.emplace_back(Bare_point(1, -1, -1), 4.);
+  points.emplace_back(Bare_point(-1, 1, -1), 4.);
+  points.emplace_back(Bare_point(-1, -1, 1), 4.);
+  points.emplace_back(Bare_point(1, 1, 1),   4.);
+  points.emplace_back(Bare_point(2, 2, 2),   1.);
+
+  // ----------------------------------------------------------------------------
+  // Init of an alpha complex from the list of points
+  // ----------------------------------------------------------------------------
+  Gudhi::alpha_complex::Alpha_complex<Kernel, true> alpha_complex_from_weighted_points(points);
+
+  Gudhi::Simplex_tree<> simplex;
+  if (alpha_complex_from_weighted_points.create_complex(simplex)) {
+    // ----------------------------------------------------------------------------
+    // Display information about the alpha complex
+    // ----------------------------------------------------------------------------
+    std::clog << "Weighted alpha complex is of dimension " << simplex.dimension() <<
+        " - " << simplex.num_simplices() << " simplices - " <<
+        simplex.num_vertices() << " vertices." << std::endl;
+
+    std::clog << "Iterator on weighted alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
+    for (auto f_simplex : simplex.filtration_simplex_range()) {
+      std::clog << "   ( ";
+      for (auto vertex : simplex.simplex_vertex_range(f_simplex)) {
+        std::clog << vertex << " ";
+      }
+      std::clog << ") -> " << "[" << simplex.filtration(f_simplex) << "] ";
+      std::clog << std::endl;
+    }
+  }
+  return 0;
+}
diff --git a/src/Alpha_complex/example/weightedalpha3dfrompoints_for_doc.txt b/src/Alpha_complex/example/weightedalpha3dfrompoints_for_doc.txt
index 7a09998d..f0695f1a 100644
--- a/src/Alpha_complex/example/weightedalpha3dfrompoints_for_doc.txt
+++ b/src/Alpha_complex/example/weightedalpha3dfrompoints_for_doc.txt
@@ -1,5 +1,5 @@
-Alpha complex is of dimension 3 - 29 simplices - 5 vertices.
-Iterator on alpha complex simplices in the filtration order, with [filtration value]:
+Weighted alpha complex is of dimension 3 - 29 simplices - 5 vertices.
+Iterator on weighted alpha complex simplices in the filtration order, with [filtration value]:
    ( 0 ) -> [-4] 
    ( 1 ) -> [-4] 
    ( 2 ) -> [-4] 
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h
index ba91998d..b315fa99 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h
@@ -12,14 +12,17 @@
 #ifndef ALPHA_COMPLEX_H_
 #define ALPHA_COMPLEX_H_
 
+#include <gudhi/Alpha_complex/Alpha_kernel_d.h>
 #include <gudhi/Debug_utils.h>
 // to construct Alpha_complex from a OFF file of points
 #include <gudhi/Points_off_io.h>
 
 #include <stdlib.h>
 #include <math.h>  // isnan, fmax
+#include <memory>  // for std::unique_ptr
 
 #include <CGAL/Delaunay_triangulation.h>
+#include <CGAL/Regular_triangulation.h>  // aka. Weighted Delaunay triangulation
 #include <CGAL/Epeck_d.h>  // For EXACT or SAFE version
 #include <CGAL/Epick_d.h>  // For FAST version
 #include <CGAL/Spatial_sort_traits_adapter_d.h>
@@ -29,6 +32,10 @@
 
 #include <Eigen/src/Core/util/Macros.h>  // for EIGEN_VERSION_AT_LEAST
 
+#include <boost/range/size.hpp>
+#include <boost/range/combine.hpp>
+#include <boost/range/adaptor/transformed.hpp>
+
 #include <iostream>
 #include <vector>
 #include <string>
@@ -91,49 +98,56 @@ template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool val
  * 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 \ref cgal &ge; 5.0.0.
  */
-template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>
+template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>, bool Weighted = false>
 class Alpha_complex {
  public:
+  /** \brief Geometric traits class that provides the geometric types and predicates needed by the triangulations.*/
+  using Geom_traits = std::conditional_t<Weighted, CGAL::Regular_triangulation_traits_adapter<Kernel>, Kernel>;
+
   // Add an int in TDS to save point index in the structure
-  typedef CGAL::Triangulation_data_structure<typename Kernel::Dimension,
-                              CGAL::Triangulation_vertex<Kernel, std::ptrdiff_t>,
-                              CGAL::Triangulation_full_cell<Kernel> > TDS;
-  /** \brief A Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/
-  typedef CGAL::Delaunay_triangulation<Kernel, TDS> Delaunay_triangulation;
-
-  /** \brief A point in Euclidean space.*/
-  typedef typename Kernel::Point_d Point_d;
-  /** \brief Geometric traits class that provides the geometric types and predicates needed by Delaunay
-   * triangulations.*/
-  typedef Kernel Geom_traits;
+  using TDS = CGAL::Triangulation_data_structure<typename Geom_traits::Dimension,
+                                                 CGAL::Triangulation_vertex<Geom_traits, std::ptrdiff_t>,
+                                                 CGAL::Triangulation_full_cell<Geom_traits> >;
 
- private:
-  typedef typename Kernel::Compute_squared_radius_d Squared_Radius;
-  typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel;
-  typedef typename Kernel::Point_dimension_d        Point_Dimension;
+  /** \brief A (Weighted or not) Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/
+  using Triangulation = std::conditional_t<Weighted, CGAL::Regular_triangulation<Kernel, TDS>,
+                                                     CGAL::Delaunay_triangulation<Kernel, TDS>>;
+
+  /** \brief CGAL kernel container for computations in function of the weighted or not characteristics.*/
+  using A_kernel_d = Alpha_kernel_d<Kernel, Weighted>;
+
+  // Numeric type of coordinates in the kernel
+  using FT = typename A_kernel_d::FT;
 
+  /** \brief Sphere is a std::pair<Kernel::Point_d, Kernel::FT> (aka. circurmcenter and squared radius).
+   * If Weighted, Sphere is a Kernel::Weighted_point_d (aka. circurmcenter and the weight value is the squared radius).
+  */
+  using Sphere = typename A_kernel_d::Sphere;
+
+  /** \brief A point, or a weighted point in Euclidean space.*/
+  using Point_d = typename Geom_traits::Point_d;
+
+ private:
   // Vertex_iterator type from CGAL.
-  typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator;
+  using CGAL_vertex_iterator = typename Triangulation::Vertex_iterator;
 
   // size_type type from CGAL.
-  typedef typename Delaunay_triangulation::size_type size_type;
+  using size_type = typename Triangulation::size_type;
 
   // Structure to switch from simplex tree vertex handle to CGAL vertex iterator.
-  typedef typename std::vector< CGAL_vertex_iterator > Vector_vertex_iterator;
-
-  // Numeric type of coordinates in the kernel
-  typedef typename Kernel::FT FT;
+  using Vector_vertex_iterator = std::vector< CGAL_vertex_iterator >;
 
  private:
   /** \brief Vertex iterator vector to switch from simplex tree vertex handle to CGAL vertex iterator.
    * Vertex handles are inserted sequentially, starting at 0.*/
   Vector_vertex_iterator vertex_handle_to_iterator_;
   /** \brief Pointer on the CGAL Delaunay triangulation.*/
-  Delaunay_triangulation* triangulation_;
+  std::unique_ptr<Triangulation> triangulation_;
   /** \brief Kernel for triangulation_ functions access.*/
-  Kernel kernel_;
+  A_kernel_d kernel_;
+
   /** \brief Cache for geometric constructions: circumcenter and squared radius of a simplex.*/
-  std::vector<std::pair<Point_d, FT>> cache_, old_cache_;
+  std::vector<Sphere> cache_, old_cache_;
 
  public:
   /** \brief Alpha_complex constructor from an OFF file name.
@@ -145,8 +159,7 @@ class Alpha_complex {
    *
    * @param[in] off_file_name OFF file [path and] name.
    */
-  Alpha_complex(const std::string& off_file_name)
-      : triangulation_(nullptr) {
+  Alpha_complex(const std::string& off_file_name) {
     Gudhi::Points_off_reader<Point_d> off_reader(off_file_name);
     if (!off_reader.is_valid()) {
       std::cerr << "Alpha_complex - Unable to read file " << off_file_name << "\n";
@@ -158,23 +171,40 @@ class Alpha_complex {
 
   /** \brief Alpha_complex constructor from a list of points.
    *
-   * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous.
+   * The vertices may be not contiguous as some points may be discarded in the triangulation (duplicate points,
+   * weighted hidden point, ...).
    * 
-   * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d
+   * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d or Kernel::Weighted_point_d.
    * 
-   * The type InputPointRange must be a range for which std::begin and
-   * std::end return input iterators on a Kernel::Point_d.
+   * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a
+   * Kernel::Point_d or Kernel::Weighted_point_d.
    */
   template<typename InputPointRange >
-  Alpha_complex(const InputPointRange& points)
-      : triangulation_(nullptr) {
+  Alpha_complex(const InputPointRange& points) {
     init_from_range(points);
   }
 
-  /** \brief Alpha_complex destructor deletes the Delaunay triangulation.
+  /** \brief Alpha_complex constructor from a list of points and weights.
+   *
+   * The vertices may be not contiguous as some points may be discarded in the triangulation (duplicate points,
+   * weighted hidden point, ...).
+   * 
+   * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d.
+   * 
+   * @param[in] weights Range of points weights. Weights must be in Kernel::FT.
+   * 
+   * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a
+   * Kernel::Point_d.
    */
-  ~Alpha_complex() {
-    delete triangulation_;
+  template <typename InputPointRange, typename WeightRange>
+  Alpha_complex(const InputPointRange& points, WeightRange weights) {
+    static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex");
+    // FIXME: this test is only valid if we have a forward range
+    GUDHI_CHECK(boost::size(weights) == boost::size(points),
+                std::invalid_argument("Points number in range different from weights range number"));
+    auto weighted_points = boost::range::combine(points, weights)
+      | boost::adaptors::transformed([](auto const&t){return Point_d(boost::get<0>(t), boost::get<1>(t));});
+    init_from_range(weighted_points);
   }
 
   // Forbid copy/move constructor/assignment operator
@@ -202,15 +232,17 @@ class Alpha_complex {
                 << std::endl;
     #endif
 
+#if CGAL_VERSION_NR < 1050101000
+    // Make compilation fail if weighted and CGAL < 5.1
+    static_assert(!Weighted, "Weighted Alpha_complex is only available for CGAL >= 5.1");
+#endif
+
     auto first = std::begin(points);
     auto last = std::end(points);
 
     if (first != last) {
-      // point_dimension function initialization
-      Point_Dimension point_dimension = kernel_.point_dimension_d_object();
-
-      // Delaunay triangulation is point dimension.
-      triangulation_ = new Delaunay_triangulation(point_dimension(*first));
+      // Delaunay triangulation init with point dimension.
+      triangulation_ = std::make_unique<Triangulation>(kernel_.get_dimension(*first));
 
       std::vector<Point_d> point_cloud(first, last);
 
@@ -218,18 +250,20 @@ class Alpha_complex {
       std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0),
                                           boost::counting_iterator<std::ptrdiff_t>(point_cloud.size()));
 
-      typedef boost::iterator_property_map<typename std::vector<Point_d>::iterator,
-                                           CGAL::Identity_property_map<std::ptrdiff_t>> Point_property_map;
-      typedef CGAL::Spatial_sort_traits_adapter_d<Kernel, Point_property_map> Search_traits_d;
+      using Point_property_map = boost::iterator_property_map<typename std::vector<Point_d>::iterator,
+                                                              CGAL::Identity_property_map<std::ptrdiff_t>>;
+      using Search_traits_d = CGAL::Spatial_sort_traits_adapter_d<Geom_traits, Point_property_map>;
 
       CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud)));
 
-      typename Delaunay_triangulation::Full_cell_handle hint;
+      typename Triangulation::Full_cell_handle hint;
       for (auto index : indices) {
-        typename Delaunay_triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint);
-        // Save index value as data to retrieve it after insertion
-        pos->data() = index;
-        hint = pos->full_cell();
+        typename Triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint);
+        if (pos != nullptr) {
+          // Save index value as data to retrieve it after insertion
+          pos->data() = index;
+          hint = pos->full_cell();
+        }
       }
       // --------------------------------------------------------------------------------------------
       // structure to retrieve CGAL points from vertex handle - one vertex handle per point.
@@ -270,9 +304,7 @@ class Alpha_complex {
       v.clear();
       for (auto vertex : cplx.simplex_vertex_range(s))
         v.push_back(get_point_(vertex));
-      Point_d c = kernel_.construct_circumcenter_d_object()(v.cbegin(), v.cend());
-      FT r = kernel_.squared_distance_d_object()(c, v[0]);
-      cache_.emplace_back(std::move(c), std::move(r));
+      cache_.emplace_back(kernel_.get_sphere(v.cbegin(), v.cend()));
     }
     return cache_[k];
   }
@@ -282,13 +314,13 @@ class Alpha_complex {
   auto radius(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) {
     auto k = cplx.key(s);
     if(k!=cplx.null_key())
-      return old_cache_[k].second;
+      return kernel_.get_squared_radius(old_cache_[k]);
     // Using a transform_range is slower, currently.
     thread_local std::vector<Point_d> v;
     v.clear();
     for (auto vertex : cplx.simplex_vertex_range(s))
       v.push_back(get_point_(vertex));
-    return kernel_.compute_squared_radius_d_object()(v.cbegin(), v.cend());
+    return kernel_.get_squared_radius(v.cbegin(), v.cend());
   }
 
  public:
@@ -322,9 +354,9 @@ class Alpha_complex {
                       bool exact = false,
                       bool default_filtration_value = false) {
     // From SimplicialComplexForAlpha type required to insert into a simplicial complex (with or without subfaces).
-    typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle;
-    typedef typename SimplicialComplexForAlpha::Simplex_handle Simplex_handle;
-    typedef std::vector<Vertex_handle> Vector_vertex;
+    using Vertex_handle = typename SimplicialComplexForAlpha::Vertex_handle;
+    using Simplex_handle = typename SimplicialComplexForAlpha::Simplex_handle;
+    using Vector_vertex = std::vector<Vertex_handle>;
 
     if (triangulation_ == nullptr) {
       std::cerr << "Alpha_complex cannot create_complex from a NULL triangulation\n";
@@ -368,6 +400,7 @@ class Alpha_complex {
     // --------------------------------------------------------------------------------------------
 
     if (!default_filtration_value) {
+      CGAL::NT_converter<FT, Filtration_value> cgal_converter;
       // --------------------------------------------------------------------------------------------
       // ### For i : d -> 0
       for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
@@ -378,14 +411,13 @@ class Alpha_complex {
             // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma)
             if (std::isnan(complex.filtration(f_simplex))) {
               Filtration_value alpha_complex_filtration = 0.0;
-              // No need to compute squared_radius on a single point - alpha is 0.0
-              if (f_simplex_dim > 0) {
+              // No need to compute squared_radius on a non-weighted single point - alpha is 0.0
+              if (Weighted || f_simplex_dim > 0) {
                 auto const& sqrad = radius(complex, f_simplex);
 #if CGAL_VERSION_NR >= 1050000000
                 if(exact) CGAL::exact(sqrad);
 #endif
-                CGAL::NT_converter<FT, Filtration_value> cv;
-                alpha_complex_filtration = cv(sqrad);
+                alpha_complex_filtration = cgal_converter(sqrad);
               }
               complex.assign_filtration(f_simplex, alpha_complex_filtration);
 #ifdef DEBUG_TRACES
@@ -393,7 +425,7 @@ class Alpha_complex {
 #endif  // DEBUG_TRACES
             }
             // No need to propagate further, unweighted points all have value 0
-            if (decr_dim > 1)
+            if (decr_dim > !Weighted)
               propagate_alpha_filtration(complex, f_simplex);
           }
         }
@@ -416,8 +448,8 @@ class Alpha_complex {
   template <typename SimplicialComplexForAlpha, typename Simplex_handle>
   void propagate_alpha_filtration(SimplicialComplexForAlpha& complex, Simplex_handle f_simplex) {
     // From SimplicialComplexForAlpha type required to assign filtration values.
-    typedef typename SimplicialComplexForAlpha::Filtration_value Filtration_value;
-    typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle;
+    using Filtration_value = typename SimplicialComplexForAlpha::Filtration_value;
+    using Vertex_handle = typename SimplicialComplexForAlpha::Vertex_handle;
 
     // ### Foreach Tau face of Sigma
     for (auto f_boundary : complex.boundary_simplex_range(f_simplex)) {
@@ -450,7 +482,7 @@ class Alpha_complex {
         while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; }
         Vertex_handle extra = *longiter;
         auto const& cache=get_cache(complex, f_boundary);
-        bool is_gab = kernel_.squared_distance_d_object()(cache.first, get_point_(extra)) >= cache.second;
+        bool is_gab = kernel_.is_gabriel(cache, get_point_(extra));
 #ifdef DEBUG_TRACES
         std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << extra << std::endl;
 #endif  // DEBUG_TRACES
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h b/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h
new file mode 100644
index 00000000..28d6d7a9
--- /dev/null
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h
@@ -0,0 +1,141 @@
+/*    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 ALPHA_COMPLEX_ALPHA_KERNEL_D_H_
+#define ALPHA_COMPLEX_ALPHA_KERNEL_D_H_
+
+#include <CGAL/version.h>  // for CGAL_VERSION_NR
+
+#include <Eigen/Core>  // for EIGEN_VERSION_AT_LEAST
+
+#include <utility>  // for std::make_pair
+
+// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
+#if CGAL_VERSION_NR < 1041101000
+# error Alpha_complex is only available for CGAL >= 4.11
+#endif
+
+#if !EIGEN_VERSION_AT_LEAST(3,1,0)
+# error Alpha_complex is only available for Eigen3 >= 3.1.0 installed with CGAL
+#endif
+
+namespace Gudhi {
+
+namespace alpha_complex {
+
+/**
+ * \class Alpha_kernel_d
+ * \brief Alpha complex kernel container.
+ * 
+ * \details
+ * The Alpha complex kernel container stores CGAL Kernel and dispatch basics computations in function of the weighted
+ * or not version of the Alpha complex.
+ */
+template < typename Kernel, bool Weighted = false >
+class Alpha_kernel_d {
+};
+
+// Unweighted Kernel_d version
+template < typename Kernel >
+class Alpha_kernel_d<Kernel, false> {
+ private:
+  // Kernel for functions access.
+  Kernel kernel_;
+ public:
+  // Fake type for compilation to succeed (cf. std::conditional in Alpha_complex.h)
+  using Weighted_point_d = void;
+  using Point_d = typename Kernel::Point_d;
+  // Numeric type of coordinates in the kernel
+  using FT = typename Kernel::FT;
+  // Sphere is a pair of point and squared radius.
+  using Sphere = typename std::pair<Point_d, FT>;
+
+  int get_dimension(const Point_d& p0) const {
+    return kernel_.point_dimension_d_object()(p0);
+  }
+
+  template<class PointIterator>
+  Sphere get_sphere(PointIterator begin, PointIterator end) const {
+    Point_d c = kernel_.construct_circumcenter_d_object()(begin, end);
+    FT r = kernel_.squared_distance_d_object()(c, *begin);
+    return std::make_pair(std::move(c), std::move(r));
+  }
+
+  template<class PointIterator>
+  FT get_squared_radius(PointIterator begin, PointIterator end) const {
+    return kernel_.compute_squared_radius_d_object()(begin, end);
+  }
+
+  FT get_squared_radius(const Sphere& sph) const {
+    return sph.second;
+  }
+
+  bool is_gabriel(const Sphere& circumcenter, const Point_d& point) {
+    return kernel_.squared_distance_d_object()(circumcenter.first, point) >= circumcenter.second;
+  }
+};
+
+// Weighted Kernel_d version
+template < typename Kernel >
+class Alpha_kernel_d<Kernel, true> {
+ private:
+  // Kernel for functions access.
+  Kernel kernel_;
+
+ public:
+  // Fake type for compilation to succeed (cf. std::conditional in Alpha_complex.h)
+  using Point_d = void;
+  using Weighted_point_d = typename Kernel::Weighted_point_d;
+  using Bare_point_d = typename Kernel::Point_d;
+  // Numeric type of coordinates in the kernel
+  using FT = typename Kernel::FT;
+  // Sphere is a weighted point (point + weight [= squared radius]).
+  using Sphere = Weighted_point_d;
+
+  int get_dimension(const Weighted_point_d& p0) const {
+    return kernel_.point_dimension_d_object()(p0.point());
+  }
+
+  template<class PointIterator>
+  Sphere get_sphere(PointIterator begin, PointIterator end) const {
+    // power_center_d_object has been renamed between CGAL 5.1 and 5.2
+#if CGAL_VERSION_NR < 1050200000
+    return kernel_.power_center_d_object()(begin, end);
+#else
+    return kernel_.construct_power_sphere_d_object()(begin, end);
+#endif
+  }
+
+  template<class PointIterator>
+  FT get_squared_radius(PointIterator begin, PointIterator end) const {
+    return kernel_.compute_squared_radius_smallest_orthogonal_sphere_d_object()(begin, end);
+  }
+
+  FT get_squared_radius(const Sphere& sph) const {
+    return sph.weight();
+  }
+
+  bool is_gabriel(const Sphere& circumcenter, const Weighted_point_d& point) {
+    // power_center_d_object has been renamed between CGAL 5.1 and 5.2
+#if CGAL_VERSION_NR < 1050200000
+    return kernel_.power_distance_d_object()(circumcenter, point) >= 0;
+#else
+    return kernel_.compute_power_product_d_object()(circumcenter, point) >= 0;
+#endif
+  }
+};
+
+}  // namespace alpha_complex
+
+namespace alphacomplex = alpha_complex;
+
+}  // namespace Gudhi
+
+#endif  // ALPHA_COMPLEX_ALPHA_KERNEL_D_H_
\ No newline at end of file
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
index 622b10ee..4e5fc933 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
@@ -154,8 +154,10 @@ class Alpha_complex_3d {
     using Kernel = CGAL::Periodic_3_regular_triangulation_traits_3<Predicates>;
   };
 
+ public:
   using Kernel = typename Kernel_3<Predicates, Weighted, Periodic>::Kernel;
 
+ private:
   using TdsVb = typename std::conditional<Periodic, CGAL::Periodic_3_triangulation_ds_vertex_base_3<>,
                                           CGAL::Triangulation_ds_vertex_base_3<>>::type;
 
diff --git a/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp
new file mode 100644
index 00000000..6da4c084
--- /dev/null
+++ b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp
@@ -0,0 +1,109 @@
+/*    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
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "alpha_kernel_d"
+#include <boost/test/unit_test.hpp>
+#include <boost/mpl/list.hpp>
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h>
+#include <CGAL/NT_converter.h>
+
+#include <iostream>
+#include <vector>
+#include <utility>  // for std::pair
+
+#include <gudhi/Alpha_complex/Alpha_kernel_d.h>
+#include <gudhi/Unitary_tests_utils.h>
+
+// Use dynamic_dimension_tag for the user to be able to set dimension
+typedef CGAL::Epeck_d< CGAL::Dynamic_dimension_tag > Exact_kernel_d;
+// Use static dimension_tag for the user not to be able to set dimension
+typedef CGAL::Epeck_d< CGAL::Dimension_tag<4> > Exact_kernel_s;
+// Use dynamic_dimension_tag for the user to be able to set dimension
+typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Inexact_kernel_d;
+// Use static dimension_tag for the user not to be able to set dimension
+typedef CGAL::Epick_d< CGAL::Dimension_tag<4> > Inexact_kernel_s;
+// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter
+
+typedef boost::mpl::list<Exact_kernel_d, Exact_kernel_s, Inexact_kernel_d, Inexact_kernel_s> list_of_kernel_variants;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_dimension, TestedKernel, list_of_kernel_variants) {
+  // Test for a point (weighted or not) in 4d, that the dimension is 4.
+
+  Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, false> kernel;
+  std::vector<double> p0 {0., 1., 2., 3.};
+  typename TestedKernel::Point_d p0_d(p0.begin(), p0.end());
+
+  std::clog << "Dimension is " << kernel.get_dimension(p0_d) << std::endl;
+  BOOST_CHECK(kernel.get_dimension(p0_d) == 4);
+
+  Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, true> w_kernel;
+  typename TestedKernel::Weighted_point_d w_p0_d(p0_d, 10.);
+
+  std::clog << "Dimension is " << w_kernel.get_dimension(w_p0_d) << std::endl;
+  BOOST_CHECK(w_kernel.get_dimension(w_p0_d) == 4);
+}
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_sphere, TestedKernel, list_of_kernel_variants) {
+  // Test with 5 points on a 3-sphere, that get_sphere returns the same center and squared radius
+  // for dD unweighted and for dD weighted with all weights at 0.
+
+  using Unweighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, false>;
+  // Sphere: (x-1)² + (y-1)² + z² + t² = 1
+  // At least 5 points for a 3-sphere
+  std::vector<double> p0 {1., 0., 0., 0.};
+  std::vector<double> p1 {0., 1., 0., 0.};
+  std::vector<double> p2 {1., 1., 1., 0.};
+  std::vector<double> p3 {1., 1., 0., 1.};
+  std::vector<double> p4 {1., 1., -1., 0.};
+
+  using Point_d = typename Unweighted_kernel::Point_d;
+  std::vector<Point_d> unw_pts;
+  unw_pts.emplace_back(p0.begin(), p0.end());
+  unw_pts.emplace_back(p1.begin(), p1.end());
+  unw_pts.emplace_back(p2.begin(), p2.end());
+  unw_pts.emplace_back(p3.begin(), p3.end());
+  unw_pts.emplace_back(p4.begin(), p4.end());
+
+  Unweighted_kernel kernel;
+  auto unw_sphere = kernel.get_sphere(unw_pts.cbegin(), unw_pts.cend());
+
+  std::clog << "Center is " << unw_sphere.first << " - squared radius is " << unw_sphere.second << std::endl;
+
+  using Weighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, true>;
+
+  using Weighted_point_d = typename Weighted_kernel::Weighted_point_d;
+  using Bare_point_d = typename Weighted_kernel::Bare_point_d;
+  std::vector<Weighted_point_d> w_pts;
+  w_pts.emplace_back(Bare_point_d(p0.begin(), p0.end()), 0.);
+  w_pts.emplace_back(Bare_point_d(p1.begin(), p1.end()), 0.);
+  w_pts.emplace_back(Bare_point_d(p2.begin(), p2.end()), 0.);
+  w_pts.emplace_back(Bare_point_d(p3.begin(), p3.end()), 0.);
+  w_pts.emplace_back(Bare_point_d(p4.begin(), p4.end()), 0.);
+
+  Weighted_kernel w_kernel;
+  auto w_sphere = w_kernel.get_sphere(w_pts.cbegin(), w_pts.cend());
+
+  std::clog << "Center is " << w_sphere.point() << " - squared radius is " << w_sphere.weight() << std::endl;
+
+  CGAL::NT_converter<typename Weighted_kernel::FT, double> cast_to_double;
+  // The results shall be the same with weights = 0.
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(unw_sphere.second), cast_to_double(w_sphere.weight()));
+  BOOST_CHECK(unw_sphere.first == w_sphere.point());
+
+  auto unw_sq_rd = kernel.get_squared_radius(unw_pts.cbegin(), unw_pts.cend());
+  std::clog << "Squared radius is " << unw_sq_rd << std::endl;
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(unw_sphere.second), cast_to_double(unw_sq_rd));
+  auto w_sq_rd = w_kernel.get_squared_radius(w_pts.cbegin(), w_pts.cend());
+  std::clog << "Squared radius is " << w_sq_rd << std::endl;
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(w_sphere.weight()), cast_to_double(w_sq_rd));
+}
diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt
index f38bd096..db5d840f 100644
--- a/src/Alpha_complex/test/CMakeLists.txt
+++ b/src/Alpha_complex/test/CMakeLists.txt
@@ -43,3 +43,20 @@ if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
   gudhi_add_boost_test(Weighted_periodic_alpha_complex_3d_test_unit)
 
 endif (NOT CGAL_VERSION VERSION_LESS 4.11.0)
+
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
+  add_executable ( Alpha_kernel_d_test_unit Alpha_kernel_d_unit_test.cpp )
+  target_link_libraries(Alpha_kernel_d_test_unit ${CGAL_LIBRARY})
+  if (TBB_FOUND)
+    target_link_libraries(Alpha_kernel_d_test_unit ${TBB_LIBRARIES})
+  endif()
+  gudhi_add_boost_test(Alpha_kernel_d_test_unit)
+  
+  add_executable ( Weighted_alpha_complex_test_unit Weighted_alpha_complex_unit_test.cpp )
+  target_link_libraries(Weighted_alpha_complex_test_unit ${CGAL_LIBRARY})
+  if (TBB_FOUND)
+    target_link_libraries(Weighted_alpha_complex_test_unit ${TBB_LIBRARIES})
+  endif()
+  gudhi_add_boost_test(Weighted_alpha_complex_test_unit)
+
+endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
\ No newline at end of file
diff --git a/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp
new file mode 100644
index 00000000..d267276c
--- /dev/null
+++ b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp
@@ -0,0 +1,229 @@
+/*    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
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "weighted_alpha_complex"
+#include <boost/test/unit_test.hpp>
+#include <boost/mpl/list.hpp>
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h>
+
+#include <cmath>  // float comparison
+#include <vector>
+#include <random>
+#include <array>
+#include <cmath> // for std::fabs
+
+#include <gudhi/Alpha_complex.h>
+#include <gudhi/Alpha_complex_3d.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Unitary_tests_utils.h>
+
+using list_of_exact_kernel_variants = boost::mpl::list<CGAL::Epeck_d< CGAL::Dynamic_dimension_tag >,
+                                                       CGAL::Epeck_d< CGAL::Dimension_tag<4> >
+                                                       > ;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(Zero_weighted_alpha_complex, Kernel, list_of_exact_kernel_variants) {
+  // Check that in exact mode for static dimension 4 the code for dD unweighted and for dD weighted with all weights
+  // 0 give exactly the same simplex tree (simplices and filtration values).
+
+  // Random points construction
+  using Point_d = typename Kernel::Point_d;
+  std::vector<Point_d> points;
+  std::uniform_real_distribution<double> rd_pts(-10., 10.);
+  std::random_device rand_dev;
+  std::mt19937 rand_engine(rand_dev());
+  for (int idx = 0; idx < 20; idx++) {
+    std::vector<double> point {rd_pts(rand_engine), rd_pts(rand_engine), rd_pts(rand_engine), rd_pts(rand_engine)};
+    points.emplace_back(point.begin(), point.end());
+  }
+  
+  // Alpha complex from points
+  Gudhi::alpha_complex::Alpha_complex<Kernel, false> alpha_complex_from_points(points);
+  Gudhi::Simplex_tree<> simplex;
+  Gudhi::Simplex_tree<>::Filtration_value infty = std::numeric_limits<Gudhi::Simplex_tree<>::Filtration_value>::infinity();
+  BOOST_CHECK(alpha_complex_from_points.create_complex(simplex, infty, true));
+  std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:"
+            << std::endl;
+  for (auto f_simplex : simplex.filtration_simplex_range()) {
+    std::clog << "   ( ";
+    for (auto vertex : simplex.simplex_vertex_range(f_simplex)) {
+      std::clog << vertex << " ";
+    }
+    std::clog << ") -> " << "[" << simplex.filtration(f_simplex) << "] " << std::endl;
+  }
+
+  // Alpha complex from zero weighted points
+  std::vector<typename Kernel::FT> weights(20, 0.);
+  Gudhi::alpha_complex::Alpha_complex<Kernel, true> alpha_complex_from_zero_weighted_points(points, weights);
+  Gudhi::Simplex_tree<> zw_simplex;
+  BOOST_CHECK(alpha_complex_from_zero_weighted_points.create_complex(zw_simplex, infty, true));
+
+  std::clog << "Iterator on zero weighted alpha complex simplices in the filtration order, with [filtration value]:"
+            << std::endl;
+  for (auto f_simplex : zw_simplex.filtration_simplex_range()) {
+    std::clog << "   ( ";
+    for (auto vertex : zw_simplex.simplex_vertex_range(f_simplex)) {
+      std::clog << vertex << " ";
+    }
+    std::clog << ") -> " << "[" << zw_simplex.filtration(f_simplex) << "] " << std::endl;
+  }
+
+  BOOST_CHECK(zw_simplex == simplex);
+}
+
+template <typename Point_d>
+bool cgal_3d_point_sort (Point_d a,Point_d b) {
+  if (a[0] != b[0])
+    return a[0] < b[0];
+  if (a[1] != b[1])
+    return a[1] < b[1];
+  return a[2] < b[2];
+}
+
+BOOST_AUTO_TEST_CASE(Weighted_alpha_complex_3d_comparison) {
+  // check that for random weighted 3d points in safe mode the 3D and dD codes give the same result with some tolerance
+
+  // Random points construction
+  using Kernel_dD = CGAL::Epeck_d< CGAL::Dimension_tag<3> >;
+  using Bare_point_d = typename Kernel_dD::Point_d;
+  using Weighted_point_d = typename Kernel_dD::Weighted_point_d;
+  std::vector<Weighted_point_d> w_points_d;
+
+  using Exact_weighted_alpha_complex_3d =
+    Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, true, false>;
+  using Bare_point_3 = typename Exact_weighted_alpha_complex_3d::Bare_point_3;
+  using Weighted_point_3 = typename Exact_weighted_alpha_complex_3d::Weighted_point_3;
+  std::vector<Weighted_point_3> w_points_3;
+
+  std::uniform_real_distribution<double> rd_pts(-10., 10.);
+  std::uniform_real_distribution<double> rd_wghts(-0.5, 0.5);
+  std::random_device rand_dev;
+  std::mt19937 rand_engine(rand_dev());
+  for (int idx = 0; idx < 20; idx++) {
+    std::vector<double> point {rd_pts(rand_engine), rd_pts(rand_engine), rd_pts(rand_engine)};
+    double weight = rd_wghts(rand_engine);
+    w_points_d.emplace_back(Bare_point_d(point.begin(), point.end()), weight);
+    w_points_3.emplace_back(Bare_point_3(point[0], point[1], point[2]), weight);
+  }
+
+  // Structures necessary for comparison
+  using Points = std::vector<std::array<double,3>>;
+  using Points_and_filtrations = std::map<Points, double>;
+  Points_and_filtrations pts_fltr_dD;
+  Points_and_filtrations pts_fltr_3d;
+
+  // Weighted alpha complex for dD version
+  Gudhi::alpha_complex::Alpha_complex<Kernel_dD, true> alpha_complex_dD_from_weighted_points(w_points_d);
+  Gudhi::Simplex_tree<> w_simplex_d;
+  BOOST_CHECK(alpha_complex_dD_from_weighted_points.create_complex(w_simplex_d));
+
+  std::clog << "Iterator on weighted alpha complex dD simplices in the filtration order, with [filtration value]:"
+            << std::endl;
+  for (auto f_simplex : w_simplex_d.filtration_simplex_range()) {
+    Points points;
+    for (auto vertex : w_simplex_d.simplex_vertex_range(f_simplex)) {
+      CGAL::NT_converter<Kernel_dD::RT, double> cgal_converter;
+      Bare_point_d pt = alpha_complex_dD_from_weighted_points.get_point(vertex).point();
+      points.push_back({cgal_converter(pt[0]), cgal_converter(pt[1]), cgal_converter(pt[2])});
+    }
+    std::clog << "   ( ";
+    std::sort (points.begin(), points.end());
+    for (auto point : points) {
+      std::clog << point[0] << " " << point[1] << " " << point[2] << " | ";
+    }
+    std::clog << ") -> " << "[" << w_simplex_d.filtration(f_simplex) << "] ";
+    std::clog << std::endl;
+    pts_fltr_dD[points] = w_simplex_d.filtration(f_simplex);
+  }
+
+  // Weighted alpha complex for 3D version
+  Exact_weighted_alpha_complex_3d alpha_complex_3D_from_weighted_points(w_points_3);
+  Gudhi::Simplex_tree<> w_simplex_3;
+  BOOST_CHECK(alpha_complex_3D_from_weighted_points.create_complex(w_simplex_3));
+
+  std::clog << "Iterator on weighted alpha complex 3D simplices in the filtration order, with [filtration value]:"
+            << std::endl;
+  for (auto f_simplex : w_simplex_3.filtration_simplex_range()) {
+    Points points;
+    for (auto vertex : w_simplex_3.simplex_vertex_range(f_simplex)) {
+      Bare_point_3 pt = alpha_complex_3D_from_weighted_points.get_point(vertex).point();
+      CGAL::NT_converter<Exact_weighted_alpha_complex_3d::Kernel::RT, double> cgal_converter;
+      points.push_back({cgal_converter(pt[0]), cgal_converter(pt[1]), cgal_converter(pt[2])});
+    }
+    std::clog << "   ( ";
+    std::sort (points.begin(), points.end());
+    for (auto point : points) {
+      std::clog << point[0] << " " << point[1] << " " << point[2] << " | ";
+    }
+    std::clog << ") -> " << "[" << w_simplex_3.filtration(f_simplex) << "] " << std::endl;
+    pts_fltr_3d[points] = w_simplex_d.filtration(f_simplex);
+  }
+
+  // Compares structures
+  auto d3_itr = pts_fltr_3d.begin();
+  auto dD_itr = pts_fltr_dD.begin();
+  for (; d3_itr != pts_fltr_3d.end() && dD_itr != pts_fltr_dD.end(); ++d3_itr) {
+    if (d3_itr->first != dD_itr->first) {
+      for(auto point : d3_itr->first)
+        std::clog << point[0] << " " << point[1] << " " << point[2] << " | ";
+      std::clog << " versus ";
+      for(auto point : dD_itr->first)
+        std::clog << point[0] << " " << point[1] << " " << point[2] << " | ";
+      std::clog << std::endl;
+      BOOST_CHECK(false);
+    }
+    // In safe mode, relative error is less than 1e-5 (can be changed with set_relative_precision_of_to_double)
+    if (std::fabs(d3_itr->second - dD_itr->second) > 1e-5 * (std::fabs(d3_itr->second) + std::fabs(dD_itr->second))) {
+      std::clog << d3_itr->second << " versus " << dD_itr->second << " diff "
+                << std::fabs(d3_itr->second - dD_itr->second) << std::endl;
+      BOOST_CHECK(false);
+    }
+    ++dD_itr;
+  }
+}
+
+using list_of_1d_kernel_variants = boost::mpl::list<CGAL::Epeck_d< CGAL::Dynamic_dimension_tag >,
+                                                    CGAL::Epeck_d< CGAL::Dimension_tag<1>>,
+                                                    CGAL::Epick_d< CGAL::Dynamic_dimension_tag >,
+                                                    CGAL::Epick_d< CGAL::Dimension_tag<1>>
+                                                    >;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(Weighted_alpha_complex_non_visible_points, Kernel, list_of_1d_kernel_variants) {
+  // check that for 2 closed weighted 1-d points, one with a high weight to hide the second one with a small weight,
+  // that the point with a small weight has the same high filtration value than the edge formed by the 2 points
+  using Point_d = typename Kernel::Point_d;
+  std::vector<Point_d> points;
+  std::vector<double> p1 {0.};
+  points.emplace_back(p1.begin(), p1.end());
+  // closed enough points
+  std::vector<double> p2 {0.1};
+  points.emplace_back(p2.begin(), p2.end());
+  std::vector<typename Kernel::FT> weights {100., 0.01};
+
+  Gudhi::alpha_complex::Alpha_complex<Kernel, true> alpha_complex(points, weights);
+  Gudhi::Simplex_tree<> stree;
+  BOOST_CHECK(alpha_complex.create_complex(stree));
+
+  std::clog << "Iterator on weighted alpha complex simplices in the filtration order, with [filtration value]:"
+            << std::endl;
+  for (auto f_simplex : stree.filtration_simplex_range()) {
+    std::clog << "   ( ";
+    for (auto vertex : stree.simplex_vertex_range(f_simplex)) {
+      std::clog << vertex << " ";
+    }
+    std::clog << ") -> " << "[" << stree.filtration(f_simplex) << "] " << std::endl;
+  }
+
+  BOOST_CHECK(stree.filtration(stree.find({0})) == -100.);
+  BOOST_CHECK(stree.filtration(stree.find({1})) == stree.filtration(stree.find({0, 1})));
+  BOOST_CHECK(stree.filtration(stree.find({1})) > 100000);
+}
\ No newline at end of file
diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt
index 1e7dc1dd..303bd0a6 100644
--- a/src/Alpha_complex/utilities/CMakeLists.txt
+++ b/src/Alpha_complex/utilities/CMakeLists.txt
@@ -1,6 +1,6 @@
 project(Alpha_complex_utilities)
 
-if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
+if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
   if (TARGET Boost::program_options)
     add_executable (alpha_complex_persistence alpha_complex_persistence.cpp)
     target_link_libraries(alpha_complex_persistence ${CGAL_LIBRARY} Boost::program_options)
@@ -28,7 +28,7 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
 
     install(TARGETS alpha_complex_persistence DESTINATION bin)
   endif()
-endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
+endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
 
 if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
   if (TARGET Boost::program_options)
diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
index e17831d9..e86b34e2 100644
--- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
@@ -17,19 +17,82 @@
 #include <gudhi/Persistent_cohomology.h>
 // to construct a simplex_tree from alpha complex
 #include <gudhi/Simplex_tree.h>
+#include <gudhi/Points_off_io.h>
 
 #include <iostream>
 #include <string>
 #include <limits>  // for numeric_limits
+#include <vector>
+#include <fstream>
 
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Filtration_value = Simplex_tree::Filtration_value;
 
 void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast,
-                     std::string &output_file_diag, Filtration_value &alpha_square_max_value,
+                     std::string &weight_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value,
                      int &coeff_field_characteristic, Filtration_value &min_persistence);
 
+template<class Point_d>
+std::vector<Point_d> read_off(const std::string &off_file_points) {
+  Gudhi::Points_off_reader<Point_d> off_reader(off_file_points);
+  if (!off_reader.is_valid()) {
+    std::cerr << "Alpha_complex - Unable to read file " << off_file_points << "\n";
+    exit(-1);  // ----- >>
+  }
+  return off_reader.get_point_cloud();
+}
+
+std::vector<double> read_weight_file(const std::string &weight_file) {
+  std::vector<double> weights;
+  // Read weights information from file
+  std::ifstream weights_ifstr(weight_file);
+  if (weights_ifstr.good()) {
+    double weight = 0.0;
+    // Attempt read the weight in a double format, return false if it fails
+    while (weights_ifstr >> weight) {
+      weights.push_back(weight);
+    }
+  } else {
+    std::cerr << "Unable to read weights file " << weight_file << std::endl;
+    exit(-1);
+  }
+  return weights;
+}
+
+template<class Kernel>
+Simplex_tree create_simplex_tree(const std::string &off_file_points, const std::string &weight_file,
+                                 bool exact_version, Filtration_value alpha_square_max_value) {
+  Simplex_tree stree;
+  auto points = read_off<typename Kernel::Point_d>(off_file_points);
+
+  if (weight_file != std::string()) {
+    std::vector<double> weights = read_weight_file(weight_file);
+    if (points.size() != weights.size()) {
+      std::cerr << "Alpha_complex - Inconsistency between number of points (" << points.size()
+                << ") and number of weights (" << weights.size() << ")" << "\n";
+      exit(-1);  // ----- >>
+    }
+    // Init of an alpha complex from an OFF file
+    Gudhi::alpha_complex::Alpha_complex<Kernel, true> alpha_complex_from_file(points, weights);
+
+    if (!alpha_complex_from_file.create_complex(stree, alpha_square_max_value, exact_version)) {
+      std::cerr << "Alpha complex simplicial complex creation failed." << std::endl;
+      exit(-1);
+    }
+  } else {
+    // Init of an alpha complex from an OFF file
+    Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(points);
+
+    if (!alpha_complex_from_file.create_complex(stree, alpha_square_max_value, exact_version)) {
+      std::cerr << "Alpha complex simplicial complex creation failed." << std::endl;
+      exit(-1);
+    }
+  }
+  return stree;
+}
+
 int main(int argc, char **argv) {
+  std::string weight_file;
   std::string off_file_points;
   std::string output_file_diag;
   bool exact_version = false;
@@ -38,48 +101,34 @@ int main(int argc, char **argv) {
   int coeff_field_characteristic;
   Filtration_value min_persistence;
 
-  program_options(argc, argv, off_file_points, exact_version, fast_version, output_file_diag, alpha_square_max_value,
-                  coeff_field_characteristic, min_persistence);
+  program_options(argc, argv, off_file_points, exact_version, fast_version, weight_file, output_file_diag,
+                  alpha_square_max_value, coeff_field_characteristic, min_persistence);
 
   if ((exact_version) && (fast_version)) {
     std::cerr << "You cannot set the exact and the fast version." << std::endl;
     exit(-1);
   }
 
-  Simplex_tree simplex;
+  Simplex_tree stree;
   if (fast_version) {
     // WARNING : CGAL::Epick_d is fast but not safe (unlike CGAL::Epeck_d)
     // (i.e. when the points are on a grid)
     using Fast_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
-
-    // Init of an alpha complex from an OFF file
-    Gudhi::alpha_complex::Alpha_complex<Fast_kernel> alpha_complex_from_file(off_file_points);
-
-    if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) {
-      std::cerr << "Fast Alpha complex simplicial complex creation failed." << std::endl;
-      exit(-1);
-    }
+    stree = create_simplex_tree<Fast_kernel>(off_file_points, weight_file, exact_version, alpha_square_max_value);
   } else {
     using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
-
-    // Init of an alpha complex from an OFF file
-    Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_points);
-
-    if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value, exact_version)) {
-      std::cerr << "Alpha complex simplicial complex creation failed." << std::endl;
-      exit(-1);
-    }
+    stree = create_simplex_tree<Kernel>(off_file_points, weight_file, exact_version, alpha_square_max_value);
   }
   // ----------------------------------------------------------------------------
   // Display information about the alpha complex
   // ----------------------------------------------------------------------------
-  std::clog << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
-            << " simplices - " << simplex.num_vertices() << " vertices." << std::endl;
+  std::clog << "Simplicial complex is of dimension " << stree.dimension() << " - " << stree.num_simplices()
+            << " simplices - " << stree.num_vertices() << " vertices." << std::endl;
 
-  std::clog << "Simplex_tree dim: " << simplex.dimension() << std::endl;
+  std::clog << "Simplex_tree dim: " << stree.dimension() << std::endl;
   // Compute the persistence diagram of the complex
   Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(
-      simplex);
+      stree);
   // initializes the coefficient field for homology
   pcoh.init_coefficients(coeff_field_characteristic);
 
@@ -98,7 +147,7 @@ int main(int argc, char **argv) {
 }
 
 void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast,
-                     std::string &output_file_diag, Filtration_value &alpha_square_max_value,
+                     std::string &weight_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value,
                      int &coeff_field_characteristic, Filtration_value &min_persistence) {
   namespace po = boost::program_options;
   po::options_description hidden("Hidden options");
@@ -111,6 +160,8 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool
       "To activate exact version of Alpha complex (default is false, not available if fast is set)")(
       "fast,f", po::bool_switch(&fast),
       "To activate fast version of Alpha complex (default is false, not available if exact is set)")(
+      "weight-file,w", po::value<std::string>(&weight_file)->default_value(std::string()),
+      "Name of file containing a point weights. Format is one weight per line:\n  W1\n  ...\n  Wn ")(
       "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
       "Name of file in which the persistence diagram is written. Default print in std::clog")(
       "max-alpha-square-value,r", po::value<Filtration_value>(&alpha_square_max_value)
diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md
index 527598a9..0d3c6027 100644
--- a/src/Alpha_complex/utilities/alphacomplex.md
+++ b/src/Alpha_complex/utilities/alphacomplex.md
@@ -46,6 +46,9 @@ for the Alpha complex construction.
 coefficient field Z/pZ for computing homology.

 * `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature

 to be recorded. Enter a negative value to see zero length intervals.

+* `-w [ --weight-file ]` is the path to the file containing the weights of the

+points (one value per line).

+Default version is not weighted.

 * `-e [ --exact ]` for the exact computation version.

 * `-f [ --fast ]` for the fast computation version.

 

@@ -58,6 +61,10 @@ to be recorded. Enter a negative value to see zero length intervals.
 N.B.:

 

 * Filtration values are alpha square values.

+* Weights values are explained on CGAL

+[dD Triangulations](https://doc.cgal.org/latest/Triangulation/index.html)

+and

+[Regular triangulation](https://doc.cgal.org/latest/Triangulation/index.html#title20) documentation.

 

 

 ## alpha_complex_3d_persistence ##