Merge tag 'tags/gudhi-release-3.2.0' into dfsg/latest

21 files changed, 399 insertions, 286 deletions

diff --git a/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp b/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp
index 99ad94b9..e7d85686 100644
--- a/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp
+++ b/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp
@@ -19,7 +19,7 @@ std::ofstream results_csv("results.csv");
 
 template <typename Kernel>
 void benchmark_points_on_torus_dD(const std::string& msg) {
-  std::cout << "+ " << msg << std::endl;
+  std::clog << "+ " << msg << std::endl;
 
   results_csv << "\"" << msg << "\";" << std::endl;
   results_csv << "\"nb_points\";"
@@ -29,7 +29,7 @@ void benchmark_points_on_torus_dD(const std::string& msg) {
 
   using K = CGAL::Epick_d<CGAL::Dimension_tag<3>>;
   for (int nb_points = 1000; nb_points <= 125000; nb_points *= 5) {
-    std::cout << "  Alpha complex dD on torus with " << nb_points << " points." << std::endl;
+    std::clog << "  Alpha complex dD on torus with " << nb_points << " points." << std::endl;
     std::vector<K::Point_d> points_on_torus = Gudhi::generate_points_on_torus_3D<K>(nb_points, 1.0, 0.5);
     std::vector<typename Kernel::Point_d> points;
 
@@ -41,26 +41,26 @@ void benchmark_points_on_torus_dD(const std::string& msg) {
     ac_create_clock.begin();
     Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_points(points);
     ac_create_clock.end();
-    std::cout << ac_create_clock;
+    std::clog << ac_create_clock;
 
     Gudhi::Simplex_tree<> complex;
     Gudhi::Clock st_create_clock("    benchmark_points_on_torus_dD - complex creation");
     st_create_clock.begin();
     alpha_complex_from_points.create_complex(complex);
     st_create_clock.end();
-    std::cout << st_create_clock;
+    std::clog << st_create_clock;
 
     results_csv << nb_points << ";" << complex.num_simplices() << ";" << ac_create_clock.num_seconds() << ";"
                 << st_create_clock.num_seconds() << ";" << std::endl;
 
-    std::cout << "    benchmark_points_on_torus_dD - nb simplices = " << complex.num_simplices() << std::endl;
+    std::clog << "    benchmark_points_on_torus_dD - nb simplices = " << complex.num_simplices() << std::endl;
   }
 }
 
 template <typename Alpha_complex_3d>
 void benchmark_points_on_torus_3D(const std::string& msg) {
   using K = CGAL::Epick_d<CGAL::Dimension_tag<3>>;
-  std::cout << "+ " << msg << std::endl;
+  std::clog << "+ " << msg << std::endl;
 
   results_csv << "\"" << msg << "\";" << std::endl;
   results_csv << "\"nb_points\";"
@@ -69,7 +69,7 @@ void benchmark_points_on_torus_3D(const std::string& msg) {
               << "\"complex_creation_time(sec.)\";" << std::endl;
 
   for (int nb_points = 1000; nb_points <= 125000; nb_points *= 5) {
-    std::cout << "  Alpha complex 3d on torus with " << nb_points << " points." << std::endl;
+    std::clog << "  Alpha complex 3d on torus with " << nb_points << " points." << std::endl;
     std::vector<K::Point_d> points_on_torus = Gudhi::generate_points_on_torus_3D<K>(nb_points, 1.0, 0.5);
     std::vector<typename Alpha_complex_3d::Point_3> points;
 
@@ -81,19 +81,19 @@ void benchmark_points_on_torus_3D(const std::string& msg) {
     ac_create_clock.begin();
     Alpha_complex_3d alpha_complex_from_points(points);
     ac_create_clock.end();
-    std::cout << ac_create_clock;
+    std::clog << ac_create_clock;
 
     Gudhi::Simplex_tree<> complex;
     Gudhi::Clock st_create_clock("    benchmark_points_on_torus_3D - complex creation");
     st_create_clock.begin();
     alpha_complex_from_points.create_complex(complex);
     st_create_clock.end();
-    std::cout << st_create_clock;
+    std::clog << st_create_clock;
 
     results_csv << nb_points << ";" << complex.num_simplices() << ";" << ac_create_clock.num_seconds() << ";"
                 << st_create_clock.num_seconds() << ";" << std::endl;
 
-    std::cout << "    benchmark_points_on_torus_3D - nb simplices = " << complex.num_simplices() << std::endl;
+    std::clog << "    benchmark_points_on_torus_3D - nb simplices = " << complex.num_simplices() << std::endl;
   }
 }
 
@@ -101,7 +101,7 @@ template <typename Weighted_alpha_complex_3d>
 void benchmark_weighted_points_on_torus_3D(const std::string& msg) {
   using K = CGAL::Epick_d<CGAL::Dimension_tag<3>>;
 
-  std::cout << "+ " << msg << std::endl;
+  std::clog << "+ " << msg << std::endl;
 
   results_csv << "\"" << msg << "\";" << std::endl;
   results_csv << "\"nb_points\";"
@@ -112,7 +112,7 @@ void benchmark_weighted_points_on_torus_3D(const std::string& msg) {
   CGAL::Random random(8);
 
   for (int nb_points = 1000; nb_points <= 125000; nb_points *= 5) {
-    std::cout << "  Alpha complex 3d on torus with " << nb_points << " points." << std::endl;
+    std::clog << "  Alpha complex 3d on torus with " << nb_points << " points." << std::endl;
     std::vector<K::Point_d> points_on_torus = Gudhi::generate_points_on_torus_3D<K>(nb_points, 1.0, 0.5);
 
     using Point = typename Weighted_alpha_complex_3d::Bare_point_3;
@@ -128,25 +128,25 @@ void benchmark_weighted_points_on_torus_3D(const std::string& msg) {
     ac_create_clock.begin();
     Weighted_alpha_complex_3d alpha_complex_from_points(points);
     ac_create_clock.end();
-    std::cout << ac_create_clock;
+    std::clog << ac_create_clock;
 
     Gudhi::Simplex_tree<> complex;
     Gudhi::Clock st_create_clock("    benchmark_weighted_points_on_torus_3D - complex creation");
     st_create_clock.begin();
     alpha_complex_from_points.create_complex(complex);
     st_create_clock.end();
-    std::cout << st_create_clock;
+    std::clog << st_create_clock;
 
     results_csv << nb_points << ";" << complex.num_simplices() << ";" << ac_create_clock.num_seconds() << ";"
                 << st_create_clock.num_seconds() << ";" << std::endl;
 
-    std::cout << "    benchmark_weighted_points_on_torus_3D - nb simplices = " << complex.num_simplices() << std::endl;
+    std::clog << "    benchmark_weighted_points_on_torus_3D - nb simplices = " << complex.num_simplices() << std::endl;
   }
 }
 
 template <typename Periodic_alpha_complex_3d>
 void benchmark_periodic_points(const std::string& msg) {
-  std::cout << "+ " << msg << std::endl;
+  std::clog << "+ " << msg << std::endl;
 
   results_csv << "\"" << msg << "\";" << std::endl;
   results_csv << "\"nb_points\";"
@@ -157,7 +157,7 @@ void benchmark_periodic_points(const std::string& msg) {
   CGAL::Random random(8);
 
   for (double nb_points = 10.; nb_points <= 40.; nb_points += 10.) {
-    std::cout << "  Periodic alpha complex 3d with " << nb_points * nb_points * nb_points << " points." << std::endl;
+    std::clog << "  Periodic alpha complex 3d with " << nb_points * nb_points * nb_points << " points." << std::endl;
     using Point = typename Periodic_alpha_complex_3d::Point_3;
     std::vector<Point> points;
 
@@ -174,25 +174,25 @@ void benchmark_periodic_points(const std::string& msg) {
     ac_create_clock.begin();
     Periodic_alpha_complex_3d alpha_complex_from_points(points, 0., 0., 0., nb_points, nb_points, nb_points);
     ac_create_clock.end();
-    std::cout << ac_create_clock;
+    std::clog << ac_create_clock;
 
     Gudhi::Simplex_tree<> complex;
     Gudhi::Clock st_create_clock("    benchmark_periodic_points - complex creation");
     st_create_clock.begin();
     alpha_complex_from_points.create_complex(complex);
     st_create_clock.end();
-    std::cout << st_create_clock;
+    std::clog << st_create_clock;
 
     results_csv << nb_points * nb_points * nb_points << ";" << complex.num_simplices() << ";"
                 << ac_create_clock.num_seconds() << ";" << st_create_clock.num_seconds() << ";" << std::endl;
 
-    std::cout << "    benchmark_periodic_points - nb simplices = " << complex.num_simplices() << std::endl;
+    std::clog << "    benchmark_periodic_points - nb simplices = " << complex.num_simplices() << std::endl;
   }
 }
 
 template <typename Weighted_periodic_alpha_complex_3d>
 void benchmark_weighted_periodic_points(const std::string& msg) {
-  std::cout << "+ " << msg << std::endl;
+  std::clog << "+ " << msg << std::endl;
 
   results_csv << "\"" << msg << "\";" << std::endl;
   results_csv << "\"nb_points\";"
@@ -203,7 +203,7 @@ void benchmark_weighted_periodic_points(const std::string& msg) {
   CGAL::Random random(8);
 
   for (double nb_points = 10.; nb_points <= 40.; nb_points += 10.) {
-    std::cout << "  Weighted periodic alpha complex 3d with " << nb_points * nb_points * nb_points << " points."
+    std::clog << "  Weighted periodic alpha complex 3d with " << nb_points * nb_points * nb_points << " points."
               << std::endl;
 
     using Point = typename Weighted_periodic_alpha_complex_3d::Bare_point_3;
@@ -224,19 +224,19 @@ void benchmark_weighted_periodic_points(const std::string& msg) {
     ac_create_clock.begin();
     Weighted_periodic_alpha_complex_3d alpha_complex_from_points(points, 0., 0., 0., nb_points, nb_points, nb_points);
     ac_create_clock.end();
-    std::cout << ac_create_clock;
+    std::clog << ac_create_clock;
 
     Gudhi::Simplex_tree<> complex;
     Gudhi::Clock st_create_clock("    benchmark_weighted_periodic_points - complex creation");
     st_create_clock.begin();
     alpha_complex_from_points.create_complex(complex);
     st_create_clock.end();
-    std::cout << st_create_clock;
+    std::clog << st_create_clock;
 
     results_csv << nb_points * nb_points * nb_points << ";" << complex.num_simplices() << ";"
                 << ac_create_clock.num_seconds() << ";" << st_create_clock.num_seconds() << ";" << std::endl;
 
-    std::cout << "    benchmark_weighted_periodic_points - nb simplices = " << complex.num_simplices() << std::endl;
+    std::clog << "    benchmark_weighted_periodic_points - nb simplices = " << complex.num_simplices() << std::endl;
   }
 }
 
diff --git a/src/Alpha_complex/concept/SimplicialComplexForAlpha.h b/src/Alpha_complex/concept/SimplicialComplexForAlpha.h
index 1c6c3b0c..c20c3201 100644
--- a/src/Alpha_complex/concept/SimplicialComplexForAlpha.h
+++ b/src/Alpha_complex/concept/SimplicialComplexForAlpha.h
@@ -72,6 +72,24 @@ struct SimplicialComplexForAlpha {
   /** \brief Return type of an insertion of a simplex
    */
   typedef unspecified Insertion_result_type;
+
+  /** \name Map interface
+   * Conceptually a `std::unordered_map<Simplex_handle,std::size_t>`.
+   *  @{ */
+  /** \brief Data stored for each simplex.
+   *
+   * Must be an integer type. */
+  typedef unspecified      Simplex_key;
+  /** \brief Returns a constant dummy number that is either negative,
+   * or at least as large as the number of simplices.  Suggested value: -1.  */
+  Simplex_key              null_key ();
+  /** \brief Returns the number stored for a simplex by `assign_key()`.
+   *
+   * If `assign_key()` has not been called, it must return `null_key()`.  */
+  Simplex_key              key      ( Simplex_handle sh );
+  /** \brief Store a number for a simplex, which can later be retrieved with `key()`.  */
+  void                     assign_key(Simplex_handle sh, Simplex_key n);
+  /** @} */
 };
 
 }  // namespace alpha_complex
diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h
index a8b1a106..60da7169 100644
--- a/src/Alpha_complex/doc/Intro_alpha_complex.h
+++ b/src/Alpha_complex/doc/Intro_alpha_complex.h
@@ -46,16 +46,17 @@ namespace alpha_complex {
  * \cite cgal:s-gkd-19b from CGAL as template parameter.
  *
  * \remark
- * - When an \f$\alpha\f$-complex is constructed with an infinite value of \f$ \alpha^2 \f$, the complex is a Delaunay
- * complex (with special filtration values).
+ * - When the simplicial complex is constructed with an infinite value of \f$ \alpha^2 \f$, the complex is a Delaunay
+ * complex with special filtration values. The Delaunay complex without filtration values is also available by passing
+ * `default_filtration_value=true` to `Alpha_complex::create_complex`.
  * - For people only interested in the topology of the \ref alpha_complex (for instance persistence),
  * \ref alpha_complex is equivalent to the \ref cech_complex and much smaller if you do not bound the radii.
  * \ref cech_complex can still make sense in higher dimension precisely because you can bound the radii.
- * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass exact=true to create_complex, the
- * filtration values are the exact ones converted to the filtration value type of the simplicial complex. This can be
- * very slow. If you pass exact=false (the default), the filtration values are only guaranteed to have a small
- * multiplicative error compared to the exact value, see <code><a class="el" target="_blank"
- * href="https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html">
+ * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass `exact=true` to
+ * `Alpha_complex::create_complex`, the filtration values are the exact ones converted to the filtration value type of
+ * the simplicial complex. This can be very slow. If you pass `exact=false` (the default), the filtration values are
+ * only guaranteed to have a small multiplicative error compared to the exact value, see <code>
+ * <a class="el" target="_blank" href="https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html">
  * CGAL::Lazy_exact_nt<NT>::set_relative_precision_of_to_double</a></code> for details. A drawback, when computing
  * persistence, is that an empty exact interval [10^12,10^12] may become a non-empty approximate interval
  * [10^12,10^12+10^6]. Using `CGAL::Epick_d` makes the computations slightly faster, and the combinatorics are still
diff --git a/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp
index 0e359a27..a2c85138 100644
--- a/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp
+++ b/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp
@@ -38,18 +38,18 @@ int main(int argc, char **argv) {
     // ----------------------------------------------------------------------------
     // Display information about the alpha complex
     // ----------------------------------------------------------------------------
-    std::cout << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
+    std::clog << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
               << " simplices - " << simplex.num_vertices() << " vertices." << std::endl;
 
-    std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
+    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::cout << "   ( ";
+      std::clog << "   ( ";
       for (auto vertex : simplex.simplex_vertex_range(f_simplex)) {
-        std::cout << vertex << " ";
+        std::clog << vertex << " ";
       }
-      std::cout << ") -> "
+      std::clog << ") -> "
                 << "[" << simplex.filtration(f_simplex) << "] ";
-      std::cout << std::endl;
+      std::clog << std::endl;
     }
   }
   return 0;
diff --git a/src/Alpha_complex/example/Alpha_complex_from_off.cpp b/src/Alpha_complex/example/Alpha_complex_from_off.cpp
index 220a66de..dba1710e 100644
--- a/src/Alpha_complex/example/Alpha_complex_from_off.cpp
+++ b/src/Alpha_complex/example/Alpha_complex_from_off.cpp
@@ -30,7 +30,7 @@ int main(int argc, char **argv) {
     ouput_file_stream.open(std::string(argv[3]));
     streambuffer = ouput_file_stream.rdbuf();
   } else {
-    streambuffer = std::cout.rdbuf();
+    streambuffer = std::clog.rdbuf();
   }
 
   Gudhi::Simplex_tree<> simplex;
diff --git a/src/Alpha_complex/example/Alpha_complex_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_from_points.cpp
index 6526ca3a..c79535bf 100644
--- a/src/Alpha_complex/example/Alpha_complex_from_points.cpp
+++ b/src/Alpha_complex/example/Alpha_complex_from_points.cpp
@@ -35,18 +35,18 @@ int main() {
     // ----------------------------------------------------------------------------
     // Display information about the alpha complex
     // ----------------------------------------------------------------------------
-    std::cout << "Alpha complex is of dimension " << simplex.dimension() <<
+    std::clog << "Alpha complex is of dimension " << simplex.dimension() <<
         " - " << simplex.num_simplices() << " simplices - " <<
         simplex.num_vertices() << " vertices." << std::endl;
 
-    std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
+    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::cout << "   ( ";
+      std::clog << "   ( ";
       for (auto vertex : simplex.simplex_vertex_range(f_simplex)) {
-        std::cout << vertex << " ";
+        std::clog << vertex << " ";
       }
-      std::cout << ") -> " << "[" << simplex.filtration(f_simplex) << "] ";
-      std::cout << std::endl;
+      std::clog << ") -> " << "[" << simplex.filtration(f_simplex) << "] ";
+      std::clog << std::endl;
     }
   }
   return 0;
diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt
index b0337934..2eecd50c 100644
--- a/src/Alpha_complex/example/CMakeLists.txt
+++ b/src/Alpha_complex/example/CMakeLists.txt
@@ -32,14 +32,18 @@ if (DIFF_PATH)
 
     add_test(Alpha_complex_example_from_off_60_diff_files ${DIFF_PATH}
         ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt)
+    set_tests_properties(Alpha_complex_example_from_off_60_diff_files PROPERTIES DEPENDS Alpha_complex_example_from_off_60) 
     add_test(Alpha_complex_example_from_off_32_diff_files ${DIFF_PATH}
         ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt)
+    set_tests_properties(Alpha_complex_example_from_off_32_diff_files PROPERTIES DEPENDS Alpha_complex_example_from_off_32) 
 
     add_test(Alpha_complex_example_fast_from_off_60_diff_files ${DIFF_PATH}
         ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt)
+    set_tests_properties(Alpha_complex_example_fast_from_off_60_diff_files PROPERTIES DEPENDS Alpha_complex_example_fast_from_off_60) 
     add_test(Alpha_complex_example_fast_from_off_32_diff_files ${DIFF_PATH}
         ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt)
-  endif()
+    set_tests_properties(Alpha_complex_example_fast_from_off_32_diff_files PROPERTIES DEPENDS Alpha_complex_example_fast_from_off_32) 
+endif()
 
   add_executable ( Alpha_complex_example_weighted_3d_from_points Weighted_alpha_complex_3d_from_points.cpp )
   target_link_libraries(Alpha_complex_example_weighted_3d_from_points  ${CGAL_LIBRARY})
diff --git a/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp b/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp
index f181005a..64728470 100644
--- a/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp
+++ b/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp
@@ -35,7 +35,7 @@ int main(int argc, char **argv) {
     ouput_file_stream.open(std::string(argv[3]));
     streambuffer = ouput_file_stream.rdbuf();
   } else {
-    streambuffer = std::cout.rdbuf();
+    streambuffer = std::clog.rdbuf();
   }
 
   Gudhi::Simplex_tree<> simplex;
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 fcf80802..c044194e 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
@@ -34,18 +34,18 @@ int main(int argc, char **argv) {
     // ----------------------------------------------------------------------------
     // Display information about the alpha complex
     // ----------------------------------------------------------------------------
-    std::cout << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
+    std::clog << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
               << " simplices - " << simplex.num_vertices() << " vertices." << std::endl;
 
-    std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
+    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::cout << "   ( ";
+      std::clog << "   ( ";
       for (auto vertex : simplex.simplex_vertex_range(f_simplex)) {
-        std::cout << vertex << " ";
+        std::clog << vertex << " ";
       }
-      std::cout << ") -> "
+      std::clog << ") -> "
                 << "[" << simplex.filtration(f_simplex) << "] ";
-      std::cout << std::endl;
+      std::clog << std::endl;
     }
   }
   return 0;
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h
index f2a05e95..ba91998d 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h
@@ -112,9 +112,6 @@ class Alpha_complex {
   typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel;
   typedef typename Kernel::Point_dimension_d        Point_Dimension;
 
-  // Type required to compute squared radius, or side of bounded sphere on a vector of points.
-  typedef typename std::vector<Point_d> Vector_of_CGAL_points;
-
   // Vertex_iterator type from CGAL.
   typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator;
 
@@ -124,6 +121,9 @@ class Alpha_complex {
   // 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;
+
  private:
   /** \brief Vertex iterator vector to switch from simplex tree vertex handle to CGAL vertex iterator.
    * Vertex handles are inserted sequentially, starting at 0.*/
@@ -132,6 +132,8 @@ class Alpha_complex {
   Delaunay_triangulation* triangulation_;
   /** \brief Kernel for triangulation_ functions access.*/
   Kernel kernel_;
+  /** \brief Cache for geometric constructions: circumcenter and squared radius of a simplex.*/
+  std::vector<std::pair<Point_d, FT>> cache_, old_cache_;
 
  public:
   /** \brief Alpha_complex constructor from an OFF file name.
@@ -237,7 +239,7 @@ class Alpha_complex {
       for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) {
         if (!triangulation_->is_infinite(*vit)) {
 #ifdef DEBUG_TRACES
-          std::cout << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl;
+          std::clog << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl;
 #endif  // DEBUG_TRACES
           vertex_handle_to_iterator_[vit->data()] = vit;
         }
@@ -246,18 +248,66 @@ class Alpha_complex {
     }
   }
 
+  /** \brief get_point_ returns the point corresponding to the vertex given as parameter.
+   * Only for internal use for faster access.
+   *
+   * @param[in] vertex Vertex handle of the point to retrieve.
+   * @return The point found.
+   */
+  const Point_d& get_point_(std::size_t vertex) const {
+    return vertex_handle_to_iterator_[vertex]->point();
+  }
+
+  /// Return a reference to the circumcenter and circumradius, writing them in the cache if necessary.
+  template<class SimplicialComplexForAlpha>
+  auto& get_cache(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) {
+    auto k = cplx.key(s);
+    if(k==cplx.null_key()){
+      k = cache_.size();
+      cplx.assign_key(s, 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));
+      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));
+    }
+    return cache_[k];
+  }
+
+  /// Return the circumradius, either from the old cache or computed, without writing to the cache.
+  template<class SimplicialComplexForAlpha>
+  auto radius(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) {
+    auto k = cplx.key(s);
+    if(k!=cplx.null_key())
+      return old_cache_[k].second;
+    // 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());
+  }
+
  public:
   /** \brief Inserts all Delaunay triangulation into the simplicial complex.
-   * It also computes the filtration values accordingly to the \ref createcomplexalgorithm
+   * It also computes the filtration values accordingly to the \ref createcomplexalgorithm if default_filtration_value
+   * is not set.
    *
    * \tparam SimplicialComplexForAlpha must meet `SimplicialComplexForAlpha` concept.
    * 
    * @param[in] complex SimplicialComplexForAlpha to be created.
    * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$, and there is very
-   * little point using anything else since it does not save time.
+   * little point using anything else since it does not save time. Useless if `default_filtration_value` is set to
+   * `true`.
    * @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not <a target="_blank"
    * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>.
-   * 
+   * @param[in] default_filtration_value Set this value to `true` if filtration values are not needed to be computed
+   * (will be set to `NaN`).
+   * Default value is `false` (which means compute the filtration values).
+   *
    * @return true if creation succeeds, false otherwise.
    * 
    * @pre Delaunay triangulation must be already constructed with dimension strictly greater than 0.
@@ -269,7 +319,8 @@ class Alpha_complex {
             typename Filtration_value = typename SimplicialComplexForAlpha::Filtration_value>
   bool create_complex(SimplicialComplexForAlpha& complex,
                       Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(),
-                      bool exact = false) {
+                      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;
@@ -296,19 +347,19 @@ class Alpha_complex {
            ++cit) {
         Vector_vertex vertexVector;
 #ifdef DEBUG_TRACES
-        std::cout << "Simplex_tree insertion ";
+        std::clog << "Simplex_tree insertion ";
 #endif  // DEBUG_TRACES
         for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
           if (*vit != nullptr) {
 #ifdef DEBUG_TRACES
-            std::cout << " " << (*vit)->data();
+            std::clog << " " << (*vit)->data();
 #endif  // DEBUG_TRACES
             // Vector of vertex construction for simplex_tree structure
             vertexVector.push_back((*vit)->data());
           }
         }
 #ifdef DEBUG_TRACES
-        std::cout << std::endl;
+        std::clog << std::endl;
 #endif  // DEBUG_TRACES
         // Insert each simplex and its subfaces in the simplex tree - filtration is NaN
         complex.insert_simplex_and_subfaces(vertexVector, std::numeric_limits<double>::quiet_NaN());
@@ -316,62 +367,48 @@ class Alpha_complex {
     }
     // --------------------------------------------------------------------------------------------
 
-    // --------------------------------------------------------------------------------------------
-    // Will be re-used many times
-    Vector_of_CGAL_points pointVector;
-    // ### For i : d -> 0
-    for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
-      // ### Foreach Sigma of dim i
-      for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) {
-        int f_simplex_dim = complex.dimension(f_simplex);
-        if (decr_dim == f_simplex_dim) {
-          pointVector.clear();
-#ifdef DEBUG_TRACES
-          std::cout << "Sigma of dim " << decr_dim << " is";
-#endif  // DEBUG_TRACES
-          for (auto vertex : complex.simplex_vertex_range(f_simplex)) {
-            pointVector.push_back(get_point(vertex));
-#ifdef DEBUG_TRACES
-            std::cout << " " << vertex;
-#endif  // DEBUG_TRACES
-          }
-#ifdef DEBUG_TRACES
-          std::cout << std::endl;
-#endif  // DEBUG_TRACES
-          // ### 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) {
-              // squared_radius function initialization
-              Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object();
-
-              CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv;
-              auto sqrad = squared_radius(pointVector.begin(), pointVector.end());
+    if (!default_filtration_value) {
+      // --------------------------------------------------------------------------------------------
+      // ### For i : d -> 0
+      for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
+        // ### Foreach Sigma of dim i
+        for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) {
+          int f_simplex_dim = complex.dimension(f_simplex);
+          if (decr_dim == f_simplex_dim) {
+            // ### 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) {
+                auto const& sqrad = radius(complex, f_simplex);
 #if CGAL_VERSION_NR >= 1050000000
-              if(exact) CGAL::exact(sqrad);
+                if(exact) CGAL::exact(sqrad);
 #endif
-              alpha_complex_filtration = cv(sqrad);
-            }
-            complex.assign_filtration(f_simplex, alpha_complex_filtration);
+                CGAL::NT_converter<FT, Filtration_value> cv;
+                alpha_complex_filtration = cv(sqrad);
+              }
+              complex.assign_filtration(f_simplex, alpha_complex_filtration);
 #ifdef DEBUG_TRACES
-            std::cout << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl;
+              std::clog << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl;
 #endif  // DEBUG_TRACES
+            }
+            // No need to propagate further, unweighted points all have value 0
+            if (decr_dim > 1)
+              propagate_alpha_filtration(complex, f_simplex);
           }
-          // No need to propagate further, unweighted points all have value 0
-          if (decr_dim > 1)
-            propagate_alpha_filtration(complex, f_simplex);
         }
+        old_cache_ = std::move(cache_);
+        cache_.clear();
       }
+      // --------------------------------------------------------------------------------------------
+  
+      // --------------------------------------------------------------------------------------------
+      // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension
+      complex.make_filtration_non_decreasing();
+      // Remove all simplices that have a filtration value greater than max_alpha_square
+      complex.prune_above_filtration(max_alpha_square);
+      // --------------------------------------------------------------------------------------------
     }
-    // --------------------------------------------------------------------------------------------
-
-    // --------------------------------------------------------------------------------------------
-    // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension
-    complex.make_filtration_non_decreasing();
-    // Remove all simplices that have a filtration value greater than max_alpha_square
-    complex.prune_above_filtration(max_alpha_square);
-    // --------------------------------------------------------------------------------------------
     return true;
   }
 
@@ -380,20 +417,18 @@ class Alpha_complex {
   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;
-#ifdef DEBUG_TRACES
     typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle;
-#endif  // DEBUG_TRACES
 
     // ### Foreach Tau face of Sigma
     for (auto f_boundary : complex.boundary_simplex_range(f_simplex)) {
 #ifdef DEBUG_TRACES
-      std::cout << " | --------------------------------------------------\n";
-      std::cout << " | Tau ";
+      std::clog << " | --------------------------------------------------\n";
+      std::clog << " | Tau ";
       for (auto vertex : complex.simplex_vertex_range(f_boundary)) {
-        std::cout << vertex << " ";
+        std::clog << vertex << " ";
       }
-      std::cout << "is a face of Sigma\n";
-      std::cout << " | isnan(complex.filtration(Tau)=" << std::isnan(complex.filtration(f_boundary)) << std::endl;
+      std::clog << "is a face of Sigma\n";
+      std::clog << " | isnan(complex.filtration(Tau)=" << std::isnan(complex.filtration(f_boundary)) << std::endl;
 #endif  // DEBUG_TRACES
       // ### If filt(Tau) is not NaN
       if (!std::isnan(complex.filtration(f_boundary))) {
@@ -402,37 +437,22 @@ class Alpha_complex {
                                                                              complex.filtration(f_simplex));
         complex.assign_filtration(f_boundary, alpha_complex_filtration);
 #ifdef DEBUG_TRACES
-        std::cout << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl;
+        std::clog << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl;
 #endif  // DEBUG_TRACES
         // ### Else
       } else {
-        // insert the Tau points in a vector for is_gabriel function
-        Vector_of_CGAL_points pointVector;
-#ifdef DEBUG_TRACES
-        Vertex_handle vertexForGabriel = Vertex_handle();
-#endif  // DEBUG_TRACES
-        for (auto vertex : complex.simplex_vertex_range(f_boundary)) {
-          pointVector.push_back(get_point(vertex));
-        }
-        // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function
-        Point_d point_for_gabriel;
-        for (auto vertex : complex.simplex_vertex_range(f_simplex)) {
-          point_for_gabriel = get_point(vertex);
-          if (std::find(pointVector.begin(), pointVector.end(), point_for_gabriel) == pointVector.end()) {
-#ifdef DEBUG_TRACES
-            // vertex is not found in Tau
-            vertexForGabriel = vertex;
-#endif  // DEBUG_TRACES
-            // No need to continue loop
-            break;
-          }
-        }
-        // is_gabriel function initialization
-        Is_Gabriel is_gabriel = kernel_.side_of_bounded_sphere_d_object();
-        bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), point_for_gabriel)
-          != CGAL::ON_BOUNDED_SIDE;
+        // Find which vertex of f_simplex is missing in f_boundary. We could actually write a variant of boundary_simplex_range that gives pairs (f_boundary, vertex). We rely on the fact that simplex_vertex_range is sorted.
+        auto longlist = complex.simplex_vertex_range(f_simplex);
+        auto shortlist = complex.simplex_vertex_range(f_boundary);
+        auto longiter = std::begin(longlist);
+        auto shortiter = std::begin(shortlist);
+        auto enditer = std::end(shortlist);
+        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;
 #ifdef DEBUG_TRACES
-        std::cout << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << vertexForGabriel << std::endl;
+        std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << extra << std::endl;
 #endif  // DEBUG_TRACES
         // ### If Tau is not Gabriel of Sigma
         if (false == is_gab) {
@@ -440,7 +460,7 @@ class Alpha_complex {
           Filtration_value alpha_complex_filtration = complex.filtration(f_simplex);
           complex.assign_filtration(f_boundary, alpha_complex_filtration);
 #ifdef DEBUG_TRACES
-          std::cout << " | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl;
+          std::clog << " | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl;
 #endif  // DEBUG_TRACES
         }
       }
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
index 7f96c94c..c19ebb79 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h
@@ -61,10 +61,7 @@ namespace Gudhi {
 
 namespace alpha_complex {
 
-#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL
-thread_local
-#endif  // GUDHI_CAN_USE_CXX11_THREAD_LOCAL
-    double RELATIVE_PRECISION_OF_TO_DOUBLE = 0.00001;
+thread_local double RELATIVE_PRECISION_OF_TO_DOUBLE = 0.00001;
 
 // Value_from_iterator returns the filtration value from an iterator on alpha shapes values
 //
@@ -472,7 +469,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
 
     alpha_shape_3_ptr_->filtration_with_alpha_values(dispatcher);
 #ifdef DEBUG_TRACES
-    std::cout << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl;
+    std::clog << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl;
 #endif  // DEBUG_TRACES
 
     using Alpha_value_iterator = typename std::vector<FT>::const_iterator;
@@ -484,7 +481,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
       if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
         for (auto i = 0; i < 4; i++) {
 #ifdef DEBUG_TRACES
-          std::cout << "from cell[" << i << "] - Point coordinates (" << (*cell)->vertex(i)->point() << ")"
+          std::clog << "from cell[" << i << "] - Point coordinates (" << (*cell)->vertex(i)->point() << ")"
                     << std::endl;
 #endif  // DEBUG_TRACES
           vertex_list.push_back((*cell)->vertex(i));
@@ -496,7 +493,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
         for (auto i = 0; i < 4; i++) {
           if ((*facet).second != i) {
 #ifdef DEBUG_TRACES
-            std::cout << "from facet=[" << i << "] - Point coordinates (" << (*facet).first->vertex(i)->point() << ")"
+            std::clog << "from facet=[" << i << "] - Point coordinates (" << (*facet).first->vertex(i)->point() << ")"
                       << std::endl;
 #endif  // DEBUG_TRACES
             vertex_list.push_back((*facet).first->vertex(i));
@@ -508,7 +505,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
       } else if (const Edge* edge = CGAL::object_cast<Edge>(&object_iterator)) {
         for (auto i : {(*edge).second, (*edge).third}) {
 #ifdef DEBUG_TRACES
-          std::cout << "from edge[" << i << "] - Point coordinates (" << (*edge).first->vertex(i)->point() << ")"
+          std::clog << "from edge[" << i << "] - Point coordinates (" << (*edge).first->vertex(i)->point() << ")"
                     << std::endl;
 #endif  // DEBUG_TRACES
           vertex_list.push_back((*edge).first->vertex(i));
@@ -519,7 +516,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
       } else if (const Alpha_vertex_handle* vertex = CGAL::object_cast<Alpha_vertex_handle>(&object_iterator)) {
 #ifdef DEBUG_TRACES
         count_vertices++;
-        std::cout << "from vertex - Point coordinates (" << (*vertex)->point() << ")" << std::endl;
+        std::clog << "from vertex - Point coordinates (" << (*vertex)->point() << ")" << std::endl;
 #endif  // DEBUG_TRACES
         vertex_list.push_back((*vertex));
       }
@@ -531,7 +528,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
           // alpha shape not found
           Complex_vertex_handle vertex = map_cgal_simplex_tree.size();
 #ifdef DEBUG_TRACES
-          std::cout << "Point (" << the_alpha_shape_vertex->point() << ") not found - insert new vertex id " << vertex
+          std::clog << "Point (" << the_alpha_shape_vertex->point() << ") not found - insert new vertex id " << vertex
                     << std::endl;
 #endif  // DEBUG_TRACES
           the_simplex.push_back(vertex);
@@ -540,7 +537,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
           // alpha shape found
           Complex_vertex_handle vertex = the_map_iterator->second;
 #ifdef DEBUG_TRACES
-          std::cout << "Point (" << the_alpha_shape_vertex->point() << ") found as vertex id " << vertex << std::endl;
+          std::clog << "Point (" << the_alpha_shape_vertex->point() << ") found as vertex id " << vertex << std::endl;
 #endif  // DEBUG_TRACES
           the_simplex.push_back(vertex);
         }
@@ -549,7 +546,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
       Filtration_value filtr = Value_from_iterator<Complexity>::perform(alpha_value_iterator);
 
 #ifdef DEBUG_TRACES
-      std::cout << "filtration = " << filtr << std::endl;
+      std::clog << "filtration = " << filtr << std::endl;
 #endif  // DEBUG_TRACES
       complex.insert_simplex(the_simplex, static_cast<Filtration_value>(filtr));
       GUDHI_CHECK(alpha_value_iterator != alpha_values.end(), "CGAL provided more simplices than values");
@@ -557,10 +554,10 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_
     }
 
 #ifdef DEBUG_TRACES
-    std::cout << "vertices \t" << count_vertices << std::endl;
-    std::cout << "edges \t\t" << count_edges << std::endl;
-    std::cout << "facets \t\t" << count_facets << std::endl;
-    std::cout << "cells \t\t" << count_cells << std::endl;
+    std::clog << "vertices \t" << count_vertices << std::endl;
+    std::clog << "edges \t\t" << count_edges << std::endl;
+    std::clog << "facets \t\t" << count_facets << std::endl;
+    std::clog << "cells \t\t" << count_cells << std::endl;
 #endif  // DEBUG_TRACES
     // --------------------------------------------------------------------------------------------
     // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension
diff --git a/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp
index cd698a27..a4ecb6ad 100644
--- a/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp
+++ b/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp
@@ -54,7 +54,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) {
   // -----------------
   // Fast version
   // -----------------
-  std::cout << "Fast alpha complex 3d" << std::endl;
+  std::clog << "Fast alpha complex 3d" << std::endl;
 
   std::vector<Fast_alpha_complex_3d::Bare_point_3> points = get_points<Fast_alpha_complex_3d::Bare_point_3>();
   Fast_alpha_complex_3d alpha_complex(points);
@@ -79,7 +79,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) {
   // -----------------
   // Exact version
   // -----------------
-  std::cout << "Exact alpha complex 3d" << std::endl;
+  std::clog << "Exact alpha complex 3d" << std::endl;
 
   std::vector<Exact_alpha_complex_3d::Bare_point_3> exact_points = get_points<Exact_alpha_complex_3d::Bare_point_3>();
   Exact_alpha_complex_3d exact_alpha_complex(exact_points);
@@ -105,13 +105,13 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) {
   // ---------------------
   // Compare both versions
   // ---------------------
-  std::cout << "Exact Alpha complex 3d is of dimension " << exact_stree.dimension() << " - Fast is "
+  std::clog << "Exact Alpha complex 3d is of dimension " << exact_stree.dimension() << " - Fast is "
             << stree.dimension() << std::endl;
   BOOST_CHECK(exact_stree.dimension() == stree.dimension());
-  std::cout << "Exact Alpha complex 3d num_simplices " << exact_stree.num_simplices() << " - Fast is "
+  std::clog << "Exact Alpha complex 3d num_simplices " << exact_stree.num_simplices() << " - Fast is "
             << stree.num_simplices() << std::endl;
   BOOST_CHECK(exact_stree.num_simplices() == stree.num_simplices());
-  std::cout << "Exact Alpha complex 3d num_vertices " << exact_stree.num_vertices() << " - Fast is "
+  std::clog << "Exact Alpha complex 3d num_vertices " << exact_stree.num_vertices() << " - Fast is "
             << stree.num_vertices() << std::endl;
   BOOST_CHECK(exact_stree.num_vertices() == stree.num_vertices());
 
@@ -119,18 +119,18 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) {
   while (sh != stree.filtration_simplex_range().end()) {
     std::vector<int> simplex;
     std::vector<int> exact_simplex;
-    std::cout << "Fast ( ";
+    std::clog << "Fast ( ";
     for (auto vertex : stree.simplex_vertex_range(*sh)) {
       simplex.push_back(vertex);
-      std::cout << vertex << " ";
+      std::clog << vertex << " ";
     }
-    std::cout << ") -> [" << stree.filtration(*sh) << "] ";
+    std::clog << ") -> [" << stree.filtration(*sh) << "] ";
 
     // Find it in the exact structure
     auto sh_exact = exact_stree.find(simplex);
     BOOST_CHECK(sh_exact != exact_stree.null_simplex());
 
-    std::cout << " versus [" << exact_stree.filtration(sh_exact) << "] " << std::endl;
+    std::clog << " versus [" << exact_stree.filtration(sh_exact) << "] " << std::endl;
     // Exact and non-exact version is not exactly the same due to float comparison
     GUDHI_TEST_FLOAT_EQUALITY_CHECK(exact_stree.filtration(sh_exact), stree.filtration(*sh));
 
@@ -139,7 +139,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) {
   // -----------------
   // Safe version
   // -----------------
-  std::cout << "Safe alpha complex 3d" << std::endl;
+  std::clog << "Safe alpha complex 3d" << std::endl;
 
   std::vector<Safe_alpha_complex_3d::Bare_point_3> safe_points = get_points<Safe_alpha_complex_3d::Bare_point_3>();
   Safe_alpha_complex_3d safe_alpha_complex(safe_points);
@@ -165,13 +165,13 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) {
   // ---------------------
   // Compare both versions
   // ---------------------
-  std::cout << "Safe Alpha complex 3d is of dimension " << safe_stree.dimension() << " - Fast is "
+  std::clog << "Safe Alpha complex 3d is of dimension " << safe_stree.dimension() << " - Fast is "
             << stree.dimension() << std::endl;
   BOOST_CHECK(safe_stree.dimension() == stree.dimension());
-  std::cout << "Safe Alpha complex 3d num_simplices " << safe_stree.num_simplices() << " - Fast is "
+  std::clog << "Safe Alpha complex 3d num_simplices " << safe_stree.num_simplices() << " - Fast is "
             << stree.num_simplices() << std::endl;
   BOOST_CHECK(safe_stree.num_simplices() == stree.num_simplices());
-  std::cout << "Safe Alpha complex 3d num_vertices " << safe_stree.num_vertices() << " - Fast is "
+  std::clog << "Safe Alpha complex 3d num_vertices " << safe_stree.num_vertices() << " - Fast is "
             << stree.num_vertices() << std::endl;
   BOOST_CHECK(safe_stree.num_vertices() == stree.num_vertices());
 
@@ -179,18 +179,18 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) {
   while (safe_sh != stree.filtration_simplex_range().end()) {
     std::vector<int> simplex;
     std::vector<int> exact_simplex;
-    std::cout << "Fast ( ";
+    std::clog << "Fast ( ";
     for (auto vertex : stree.simplex_vertex_range(*safe_sh)) {
       simplex.push_back(vertex);
-      std::cout << vertex << " ";
+      std::clog << vertex << " ";
     }
-    std::cout << ") -> [" << stree.filtration(*safe_sh) << "] ";
+    std::clog << ") -> [" << stree.filtration(*safe_sh) << "] ";
 
     // Find it in the exact structure
     auto sh_exact = safe_stree.find(simplex);
     BOOST_CHECK(sh_exact != safe_stree.null_simplex());
 
-    std::cout << " versus [" << safe_stree.filtration(sh_exact) << "] " << std::endl;
+    std::clog << " versus [" << safe_stree.filtration(sh_exact) << "] " << std::endl;
     // Exact and non-exact version is not exactly the same due to float comparison
     GUDHI_TEST_FLOAT_EQUALITY_CHECK(safe_stree.filtration(sh_exact), stree.filtration(*safe_sh), 1e-15);
 
diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
index 27b671dd..4b37e4bd 100644
--- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
+++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
@@ -48,7 +48,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of
   // ----------------------------------------------------------------------------
   std::string off_file_name("alphacomplexdoc.off");
   double max_alpha_square_value = 60.0;
-  std::cout << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" <<
+  std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" <<
       max_alpha_square_value << "==========" << std::endl;
 
   Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_file(off_file_name);
@@ -56,29 +56,29 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of
   Gudhi::Simplex_tree<> simplex_tree_60;
   BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_60, max_alpha_square_value));
 
-  std::cout << "simplex_tree_60.dimension()=" << simplex_tree_60.dimension() << std::endl;
+  std::clog << "simplex_tree_60.dimension()=" << simplex_tree_60.dimension() << std::endl;
   BOOST_CHECK(simplex_tree_60.dimension() == 2);
 
-  std::cout << "simplex_tree_60.num_vertices()=" << simplex_tree_60.num_vertices() << std::endl;
+  std::clog << "simplex_tree_60.num_vertices()=" << simplex_tree_60.num_vertices() << std::endl;
   BOOST_CHECK(simplex_tree_60.num_vertices() == 7);
 
-  std::cout << "simplex_tree_60.num_simplices()=" << simplex_tree_60.num_simplices() << std::endl;
+  std::clog << "simplex_tree_60.num_simplices()=" << simplex_tree_60.num_simplices() << std::endl;
   BOOST_CHECK(simplex_tree_60.num_simplices() == 25);
 
   max_alpha_square_value = 59.0;
-  std::cout << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" <<
+  std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" <<
       max_alpha_square_value << "==========" << std::endl;
 
   Gudhi::Simplex_tree<> simplex_tree_59;
   BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_59, max_alpha_square_value));
   
-  std::cout << "simplex_tree_59.dimension()=" << simplex_tree_59.dimension() << std::endl;
+  std::clog << "simplex_tree_59.dimension()=" << simplex_tree_59.dimension() << std::endl;
   BOOST_CHECK(simplex_tree_59.dimension() == 2);
 
-  std::cout << "simplex_tree_59.num_vertices()=" << simplex_tree_59.num_vertices() << std::endl;
+  std::clog << "simplex_tree_59.num_vertices()=" << simplex_tree_59.num_vertices() << std::endl;
   BOOST_CHECK(simplex_tree_59.num_vertices() == 7);
 
-  std::cout << "simplex_tree_59.num_simplices()=" << simplex_tree_59.num_simplices() << std::endl;
+  std::clog << "simplex_tree_59.num_simplices()=" << simplex_tree_59.num_simplices() << std::endl;
   BOOST_CHECK(simplex_tree_59.num_simplices() == 23);
 }
 
@@ -115,30 +115,30 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
   // ----------------------------------------------------------------------------
   Gudhi::alpha_complex::Alpha_complex<Kernel_4> alpha_complex_from_points(points);
 
-  std::cout << "========== Alpha_complex_from_points ==========" << std::endl;
+  std::clog << "========== Alpha_complex_from_points ==========" << std::endl;
 
   Gudhi::Simplex_tree<> simplex_tree;
   BOOST_CHECK(alpha_complex_from_points.create_complex(simplex_tree));
   
   // Another way to check num_simplices
-  std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
+  std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
   int num_simplices = 0;
   for (auto f_simplex : simplex_tree.filtration_simplex_range()) {
     num_simplices++;
-    std::cout << "   ( ";
+    std::clog << "   ( ";
     for (auto vertex : simplex_tree.simplex_vertex_range(f_simplex)) {
-      std::cout << vertex << " ";
+      std::clog << vertex << " ";
     }
-    std::cout << ") -> " << "[" << simplex_tree.filtration(f_simplex) << "] ";
-    std::cout << std::endl;
+    std::clog << ") -> " << "[" << simplex_tree.filtration(f_simplex) << "] ";
+    std::clog << std::endl;
   }
   BOOST_CHECK(num_simplices == 15);
-  std::cout << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl;
+  std::clog << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl;
   BOOST_CHECK(simplex_tree.num_simplices() == 15);
 
-  std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl;
+  std::clog << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl;
   BOOST_CHECK(simplex_tree.dimension() == 3);
-  std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl;
+  std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl;
   BOOST_CHECK(simplex_tree.num_vertices() == points.size());
 
   for (auto f_simplex : simplex_tree.filtration_simplex_range()) {
@@ -162,22 +162,22 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
   }
 
   Point_4 p0 = alpha_complex_from_points.get_point(0);
-  std::cout << "alpha_complex_from_points.get_point(0)=" << p0 << std::endl;
+  std::clog << "alpha_complex_from_points.get_point(0)=" << p0 << std::endl;
   BOOST_CHECK(4 == p0.dimension());
   BOOST_CHECK(is_point_in_list(points, p0));
 
   Point_4 p1 = alpha_complex_from_points.get_point(1);
-  std::cout << "alpha_complex_from_points.get_point(1)=" << p1 << std::endl;
+  std::clog << "alpha_complex_from_points.get_point(1)=" << p1 << std::endl;
   BOOST_CHECK(4 == p1.dimension());
   BOOST_CHECK(is_point_in_list(points, p1));
 
   Point_4 p2 = alpha_complex_from_points.get_point(2);
-  std::cout << "alpha_complex_from_points.get_point(2)=" << p2 << std::endl;
+  std::clog << "alpha_complex_from_points.get_point(2)=" << p2 << std::endl;
   BOOST_CHECK(4 == p2.dimension());
   BOOST_CHECK(is_point_in_list(points, p2));
 
   Point_4 p3 = alpha_complex_from_points.get_point(3);
-  std::cout << "alpha_complex_from_points.get_point(3)=" << p3 << std::endl;
+  std::clog << "alpha_complex_from_points.get_point(3)=" << p3 << std::endl;
   BOOST_CHECK(4 == p3.dimension());
   BOOST_CHECK(is_point_in_list(points, p3));
 
@@ -188,30 +188,27 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
   
   // Test after prune_above_filtration
   bool modified = simplex_tree.prune_above_filtration(0.6);
-  if (modified) {
-    simplex_tree.initialize_filtration();
-  }
   BOOST_CHECK(modified);
   
   // Another way to check num_simplices
-  std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
+  std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl;
   num_simplices = 0;
   for (auto f_simplex : simplex_tree.filtration_simplex_range()) {
     num_simplices++;
-    std::cout << "   ( ";
+    std::clog << "   ( ";
     for (auto vertex : simplex_tree.simplex_vertex_range(f_simplex)) {
-      std::cout << vertex << " ";
+      std::clog << vertex << " ";
     }
-    std::cout << ") -> " << "[" << simplex_tree.filtration(f_simplex) << "] ";
-    std::cout << std::endl;
+    std::clog << ") -> " << "[" << simplex_tree.filtration(f_simplex) << "] ";
+    std::clog << std::endl;
   }
   BOOST_CHECK(num_simplices == 10);
-  std::cout << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl;
+  std::clog << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl;
   BOOST_CHECK(simplex_tree.num_simplices() == 10);
 
-  std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl;
+  std::clog << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl;
   BOOST_CHECK(simplex_tree.dimension() == 1);
-  std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl;
+  std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl;
   BOOST_CHECK(simplex_tree.num_vertices() == 4);
 
   for (auto f_simplex : simplex_tree.filtration_simplex_range()) {
@@ -231,7 +228,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
 }
 
 BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, list_of_kernel_variants) {
-  std::cout << "========== Alpha_complex_from_empty_points ==========" << std::endl;
+  std::clog << "========== Alpha_complex_from_empty_points ==========" << std::endl;
 
   // ----------------------------------------------------------------------------
   // Init of an empty list of points
@@ -249,13 +246,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, lis
   Gudhi::Simplex_tree<> simplex_tree;
   BOOST_CHECK(!alpha_complex_from_points.create_complex(simplex_tree));
   
-  std::cout << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl;
+  std::clog << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl;
   BOOST_CHECK(simplex_tree.num_simplices() == 0);
 
-  std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl;
+  std::clog << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl;
   BOOST_CHECK(simplex_tree.dimension() == -1);
   
-  std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl;
+  std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl;
   BOOST_CHECK(simplex_tree.num_vertices() == points.size());
 }
 
@@ -264,7 +261,7 @@ using Exact_kernel_2 = CGAL::Epeck_d< CGAL::Dimension_tag<2> >;
 using list_of_kernel_2_variants = boost::mpl::list<Inexact_kernel_2, Exact_kernel_2>;
 
 BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_with_duplicated_points, TestedKernel, list_of_kernel_2_variants) {
-  std::cout << "========== Alpha_complex_with_duplicated_points ==========" << std::endl;
+  std::clog << "========== Alpha_complex_with_duplicated_points ==========" << std::endl;
 
   using Point = typename TestedKernel::Point_d;
   using Vector_of_points = std::vector<Point>;
@@ -287,14 +284,14 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_with_duplicated_points, TestedKernel
   // ----------------------------------------------------------------------------
   // Init of an alpha complex from the list of points
   // ----------------------------------------------------------------------------
-  std::cout << "Init" << std::endl;
+  std::clog << "Init" << std::endl;
   Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_points(points);
 
   Gudhi::Simplex_tree<> simplex_tree;
-  std::cout << "create_complex" << std::endl;
+  std::clog << "create_complex" << std::endl;
   BOOST_CHECK(alpha_complex_from_points.create_complex(simplex_tree));
   
-  std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices()
+  std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices()
       << std::endl;
   BOOST_CHECK(simplex_tree.num_vertices() < points.size());
 }
diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt
index 0476c6d4..fe4b23e4 100644
--- a/src/Alpha_complex/test/CMakeLists.txt
+++ b/src/Alpha_complex/test/CMakeLists.txt
@@ -8,11 +8,15 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
 
   add_executable ( Alpha_complex_test_unit Alpha_complex_unit_test.cpp )
   target_link_libraries(Alpha_complex_test_unit ${CGAL_LIBRARY})
+  add_executable ( Delaunay_complex_test_unit Delaunay_complex_unit_test.cpp )
+  target_link_libraries(Delaunay_complex_test_unit ${CGAL_LIBRARY})
   if (TBB_FOUND)
     target_link_libraries(Alpha_complex_test_unit ${TBB_LIBRARIES})
+    target_link_libraries(Delaunay_complex_test_unit ${TBB_LIBRARIES})
   endif()
 
   gudhi_add_boost_test(Alpha_complex_test_unit)
+  gudhi_add_boost_test(Delaunay_complex_test_unit)
 
   add_executable ( Alpha_complex_3d_test_unit Alpha_complex_3d_unit_test.cpp )
   target_link_libraries(Alpha_complex_3d_test_unit ${CGAL_LIBRARY})
diff --git a/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp b/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp
new file mode 100644
index 00000000..c1cc1fab
--- /dev/null
+++ b/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp
@@ -0,0 +1,68 @@
+/*    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 "delaunay_complex"
+#include <boost/test/unit_test.hpp>
+#include <boost/mpl/list.hpp>
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h>
+
+#include <vector>
+#include <limits>  // NaN
+#include <cmath>
+
+#include <gudhi/Alpha_complex.h>
+// to construct a simplex_tree from Delaunay_triangulation
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Unitary_tests_utils.h>
+#include <gudhi/random_point_generators.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<5> > 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<5> > 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;
+
+using Simplex_tree = Gudhi::Simplex_tree<>;
+using Simplex_handle = Simplex_tree::Simplex_handle;
+
+BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of_kernel_variants) {
+  std::cout << "*****************************************************************************************************";
+  using Point = typename TestedKernel::Point_d;
+  std::vector<Point> points;
+  // 50 points on a 4-sphere
+  points = Gudhi::generate_points_on_sphere_d<TestedKernel>(10, 5, 1.);
+
+  Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex(points);
+
+  // Alpha complex
+  Simplex_tree stree_from_alpha_complex;
+  BOOST_CHECK(alpha_complex.create_complex(stree_from_alpha_complex));
+
+  // Delaunay complex
+  Simplex_tree stree_from_delaunay_complex;
+  BOOST_CHECK(alpha_complex.create_complex(stree_from_delaunay_complex, 0., false, true));
+
+  // Check all the simplices from alpha complex are in the Delaunay complex
+  for (auto f_simplex : stree_from_alpha_complex.complex_simplex_range()) {
+    Simplex_handle sh = stree_from_delaunay_complex.find(stree_from_alpha_complex.simplex_vertex_range(f_simplex));
+    BOOST_CHECK(std::isnan(stree_from_delaunay_complex.filtration(sh)));
+    BOOST_CHECK(sh != stree_from_delaunay_complex.null_simplex());
+  }
+}
diff --git a/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp
index 731763fa..9eef920b 100644
--- a/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp
+++ b/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp
@@ -43,14 +43,14 @@ typedef boost::mpl::list<Fast_periodic_alpha_complex_3d, Safe_periodic_alpha_com
     periodic_variants_type_list;
 
 BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_periodic_throw, Periodic_alpha_complex_3d, periodic_variants_type_list) {
-  std::cout << "Periodic alpha complex 3d exception throw" << std::endl;
+  std::clog << "Periodic alpha complex 3d exception throw" << std::endl;
   using Bare_point_3 = typename Periodic_alpha_complex_3d::Bare_point_3;
   std::vector<Bare_point_3> p_points;
 
   // Not important, this is not what we want to check
   p_points.push_back(Bare_point_3(0.0, 0.0, 0.0));
 
-  std::cout << "Check exception throw in debug mode" << std::endl;
+  std::clog << "Check exception throw in debug mode" << std::endl;
   // Check it throws an exception when the cuboid is not iso
   BOOST_CHECK_THROW(Periodic_alpha_complex_3d periodic_alpha_complex(p_points, 0., 0., 0., 0.9, 1., 1.),
                     std::invalid_argument);
@@ -71,7 +71,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) {
   // ---------------------
   // Fast periodic version
   // ---------------------
-  std::cout << "Fast periodic alpha complex 3d" << std::endl;
+  std::clog << "Fast periodic alpha complex 3d" << std::endl;
 
   using Creator = CGAL::Creator_uniform_3<double, Fast_periodic_alpha_complex_3d::Bare_point_3>;
   CGAL::Random random(7);
@@ -106,7 +106,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) {
   // ----------------------
   // Exact periodic version
   // ----------------------
-  std::cout << "Exact periodic alpha complex 3d" << std::endl;
+  std::clog << "Exact periodic alpha complex 3d" << std::endl;
 
   std::vector<Exact_periodic_alpha_complex_3d::Bare_point_3> e_p_points;
 
@@ -122,13 +122,13 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) {
   // ---------------------
   // Compare both versions
   // ---------------------
-  std::cout << "Exact periodic alpha complex 3d is of dimension " << exact_stree.dimension() << " - Non exact is "
+  std::clog << "Exact periodic alpha complex 3d is of dimension " << exact_stree.dimension() << " - Non exact is "
             << stree.dimension() << std::endl;
   BOOST_CHECK(exact_stree.dimension() == stree.dimension());
-  std::cout << "Exact periodic alpha complex 3d num_simplices " << exact_stree.num_simplices() << " - Non exact is "
+  std::clog << "Exact periodic alpha complex 3d num_simplices " << exact_stree.num_simplices() << " - Non exact is "
             << stree.num_simplices() << std::endl;
   BOOST_CHECK(exact_stree.num_simplices() == stree.num_simplices());
-  std::cout << "Exact periodic alpha complex 3d num_vertices " << exact_stree.num_vertices() << " - Non exact is "
+  std::clog << "Exact periodic alpha complex 3d num_vertices " << exact_stree.num_vertices() << " - Non exact is "
             << stree.num_vertices() << std::endl;
   BOOST_CHECK(exact_stree.num_vertices() == stree.num_vertices());
 
@@ -155,7 +155,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) {
   // ----------------------
   // Safe periodic version
   // ----------------------
-  std::cout << "Safe periodic alpha complex 3d" << std::endl;
+  std::clog << "Safe periodic alpha complex 3d" << std::endl;
 
   std::vector<Safe_periodic_alpha_complex_3d::Bare_point_3> s_p_points;
 
diff --git a/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp
index 8035f6e8..6b31bea6 100644
--- a/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp
+++ b/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp
@@ -55,13 +55,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_throw, Weighted_alpha_compl
   // weights size is different from w_points size to make weighted Alpha_complex_3d throw in debug mode
   std::vector<double> weights = {0.01, 0.005, 0.006, 0.01, 0.009, 0.001};
 
-  std::cout << "Check exception throw in debug mode" << std::endl;
+  std::clog << "Check exception throw in debug mode" << std::endl;
   BOOST_CHECK_THROW(Weighted_alpha_complex_3d wac(w_points, weights), std::invalid_argument);
 }
 #endif
 
 BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted, Weighted_alpha_complex_3d, weighted_variants_type_list) {
-  std::cout << "Weighted alpha complex 3d from points and weights" << std::endl;
+  std::clog << "Weighted alpha complex 3d from points and weights" << std::endl;
   using Bare_point_3 = typename Weighted_alpha_complex_3d::Bare_point_3;
   std::vector<Bare_point_3> w_points;
   w_points.push_back(Bare_point_3(0.0, 0.0, 0.0));
@@ -78,7 +78,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted, Weighted_alpha_complex_3d,
   Gudhi::Simplex_tree<> stree;
   alpha_complex_p_a_w.create_complex(stree);
 
-  std::cout << "Weighted alpha complex 3d from weighted points" << std::endl;
+  std::clog << "Weighted alpha complex 3d from weighted points" << std::endl;
   using Weighted_point_3 = typename Weighted_alpha_complex_3d::Weighted_point_3;
 
   std::vector<Weighted_point_3> weighted_points;
@@ -112,13 +112,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted, Weighted_alpha_complex_3d,
   // ---------------------
   // Compare both versions
   // ---------------------
-  std::cout << "Weighted alpha complex 3d is of dimension " << stree_bis.dimension() << " - versus "
+  std::clog << "Weighted alpha complex 3d is of dimension " << stree_bis.dimension() << " - versus "
             << stree.dimension() << std::endl;
   BOOST_CHECK(stree_bis.dimension() == stree.dimension());
-  std::cout << "Weighted alpha complex 3d num_simplices " << stree_bis.num_simplices() << " - versus "
+  std::clog << "Weighted alpha complex 3d num_simplices " << stree_bis.num_simplices() << " - versus "
             << stree.num_simplices() << std::endl;
   BOOST_CHECK(stree_bis.num_simplices() == stree.num_simplices());
-  std::cout << "Weighted alpha complex 3d num_vertices " << stree_bis.num_vertices() << " - versus "
+  std::clog << "Weighted alpha complex 3d num_vertices " << stree_bis.num_vertices() << " - versus "
             << stree.num_vertices() << std::endl;
   BOOST_CHECK(stree_bis.num_vertices() == stree.num_vertices());
 
@@ -127,18 +127,18 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted, Weighted_alpha_complex_3d,
     std::vector<int> simplex;
     std::vector<int> exact_simplex;
 #ifdef DEBUG_TRACES
-    std::cout << " ( ";
+    std::clog << " ( ";
 #endif
     for (auto vertex : stree.simplex_vertex_range(*sh)) {
       simplex.push_back(vertex);
 #ifdef DEBUG_TRACES
-      std::cout << vertex << " ";
+      std::clog << vertex << " ";
 #endif
     }
 #ifdef DEBUG_TRACES
-    std::cout << ") -> "
+    std::clog << ") -> "
               << "[" << stree.filtration(*sh) << "] ";
-    std::cout << std::endl;
+    std::clog << std::endl;
 #endif
 
     // Find it in the exact structure
diff --git a/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp
index b09e92d5..610b9f3d 100644
--- a/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp
+++ b/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp
@@ -45,7 +45,7 @@ typedef boost::mpl::list<Fast_weighted_periodic_alpha_complex_3d, Exact_weighted
 #ifdef GUDHI_DEBUG
 BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_periodic_throw, Weighted_periodic_alpha_complex_3d,
                               wp_variants_type_list) {
-  std::cout << "Weighted periodic alpha complex 3d exception throw" << std::endl;
+  std::clog << "Weighted periodic alpha complex 3d exception throw" << std::endl;
 
   using Creator = CGAL::Creator_uniform_3<double, typename Weighted_periodic_alpha_complex_3d::Bare_point_3>;
   CGAL::Random random(7);
@@ -62,7 +62,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_periodic_throw, Weighted_pe
     p_weights.push_back(random.get_double(0., 0.01));
   }
 
-  std::cout << "Cuboid is not iso exception" << std::endl;
+  std::clog << "Cuboid is not iso exception" << std::endl;
   // Check it throws an exception when the cuboid is not iso
   BOOST_CHECK_THROW(
       Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, -1., -1., -1., 0.9, 1., 1.),
@@ -83,7 +83,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_periodic_throw, Weighted_pe
       Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, -1., -1., -1., 1., 1., 1.1),
       std::invalid_argument);
 
-  std::cout << "0 <= point.weight() < 1/64 * domain_size * domain_size exception" << std::endl;
+  std::clog << "0 <= point.weight() < 1/64 * domain_size * domain_size exception" << std::endl;
   // Weights must be in range ]0, 1/64 = 0.015625[
   double temp = p_weights[25];
   p_weights[25] = 1.0;
@@ -97,7 +97,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_periodic_throw, Weighted_pe
                     std::invalid_argument);
   p_weights[14] = temp;
 
-  std::cout << "wp_points and p_weights size exception" << std::endl;
+  std::clog << "wp_points and p_weights size exception" << std::endl;
   // Weights and points must have the same size
   // + 1
   p_weights.push_back(1e-10);
@@ -115,7 +115,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) {
   // ---------------------
   // Fast weighted periodic version
   // ---------------------
-  std::cout << "Fast weighted periodic alpha complex 3d" << std::endl;
+  std::clog << "Fast weighted periodic alpha complex 3d" << std::endl;
 
   using Creator = CGAL::Creator_uniform_3<double, Fast_weighted_periodic_alpha_complex_3d::Bare_point_3>;
   CGAL::Random random(7);
@@ -140,7 +140,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) {
   // ----------------------
   // Exact weighted periodic version
   // ----------------------
-  std::cout << "Exact weighted periodic alpha complex 3d" << std::endl;
+  std::clog << "Exact weighted periodic alpha complex 3d" << std::endl;
 
   std::vector<Exact_weighted_periodic_alpha_complex_3d::Bare_point_3> e_p_points;
 
@@ -156,13 +156,13 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) {
   // ---------------------
   // Compare both versions
   // ---------------------
-  std::cout << "Exact weighted periodic alpha complex 3d is of dimension " << exact_stree.dimension()
+  std::clog << "Exact weighted periodic alpha complex 3d is of dimension " << exact_stree.dimension()
             << " - Non exact is " << stree.dimension() << std::endl;
   BOOST_CHECK(exact_stree.dimension() == stree.dimension());
-  std::cout << "Exact weighted periodic alpha complex 3d num_simplices " << exact_stree.num_simplices()
+  std::clog << "Exact weighted periodic alpha complex 3d num_simplices " << exact_stree.num_simplices()
             << " - Non exact is " << stree.num_simplices() << std::endl;
   BOOST_CHECK(exact_stree.num_simplices() == stree.num_simplices());
-  std::cout << "Exact weighted periodic alpha complex 3d num_vertices " << exact_stree.num_vertices()
+  std::clog << "Exact weighted periodic alpha complex 3d num_vertices " << exact_stree.num_vertices()
             << " - Non exact is " << stree.num_vertices() << std::endl;
   BOOST_CHECK(exact_stree.num_vertices() == stree.num_vertices());
 
@@ -189,7 +189,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) {
   // ----------------------
   // Safe weighted periodic version
   // ----------------------
-  std::cout << "Safe weighted periodic alpha complex 3d" << std::endl;
+  std::clog << "Safe weighted periodic alpha complex 3d" << std::endl;
 
   std::vector<Safe_weighted_periodic_alpha_complex_3d::Bare_point_3> s_p_points;
 
diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt
index 57b92942..2ffbdde0 100644
--- a/src/Alpha_complex/utilities/CMakeLists.txt
+++ b/src/Alpha_complex/utilities/CMakeLists.txt
@@ -2,7 +2,7 @@ project(Alpha_complex_utilities)
 
 if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
   add_executable (alpha_complex_persistence alpha_complex_persistence.cpp)
-  target_link_libraries(alpha_complex_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+  target_link_libraries(alpha_complex_persistence ${CGAL_LIBRARY} Boost::program_options)
 
   if (TBB_FOUND)
     target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES})
@@ -16,14 +16,20 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
   if (DIFF_PATH)
     add_test(Alpha_complex_utilities_diff_exact_alpha_complex ${DIFF_PATH}
              "exact.pers" "safe.pers")
+    set_tests_properties(Alpha_complex_utilities_diff_exact_alpha_complex PROPERTIES DEPENDS
+             "Alpha_complex_utilities_exact_alpha_complex_persistence;Alpha_complex_utilities_safe_alpha_complex_persistence")
+             
     add_test(Alpha_complex_utilities_diff_fast_alpha_complex ${DIFF_PATH}
              "fast.pers" "safe.pers")
+    set_tests_properties(Alpha_complex_utilities_diff_fast_alpha_complex PROPERTIES DEPENDS
+             "Alpha_complex_utilities_fast_alpha_complex_persistence;Alpha_complex_utilities_safe_alpha_complex_persistence")
+
   endif()
   
   install(TARGETS alpha_complex_persistence DESTINATION bin)
 
   add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp)
-  target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
+  target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY} Boost::program_options)
   if (TBB_FOUND)
     target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES})
   endif(TBB_FOUND)
@@ -36,15 +42,19 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
           "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off"
           "-p" "2" "-m" "0.45" "-o" "exact_3d.pers" "-e")
 
-  add_test(NAME Alpha_complex_utilities_safe_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
+  add_test(NAME Alpha_complex_utilities_fast_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
           "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off"
           "-p" "2" "-m" "0.45" "-o" "fast_3d.pers" "-f")
 
   if (DIFF_PATH)
     add_test(Alpha_complex_utilities_diff_exact_alpha_complex_3d ${DIFF_PATH}
             "exact_3d.pers" "safe_3d.pers")
+    set_tests_properties(Alpha_complex_utilities_diff_exact_alpha_complex_3d PROPERTIES DEPENDS
+            "Alpha_complex_utilities_exact_alpha_complex_3d;Alpha_complex_utilities_alpha_complex_3d")
     add_test(Alpha_complex_utilities_diff_fast_alpha_complex_3d ${DIFF_PATH}
             "fast_3d.pers" "safe_3d.pers")
+    set_tests_properties(Alpha_complex_utilities_diff_fast_alpha_complex_3d PROPERTIES DEPENDS
+            "Alpha_complex_utilities_fast_alpha_complex_3d;Alpha_complex_utilities_alpha_complex_3d")
   endif()
 
   add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
index 929fc2e8..91899040 100644
--- a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
@@ -222,10 +222,7 @@ int main(int argc, char **argv) {
       break;
   }
 
-  // Sort the simplices in the order of the filtration
-  simplex_tree.initialize_filtration();
-
-  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
+  std::clog << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
   // Compute the persistence diagram of the complex
   Persistent_cohomology pcoh(simplex_tree, true);
   // initializes the coefficient field for homology
@@ -237,7 +234,7 @@ int main(int argc, char **argv) {
   if (output_file_diag.empty()) {
     pcoh.output_diagram();
   } else {
-    std::cout << "Result in file: " << output_file_diag << std::endl;
+    std::clog << "Result in file: " << output_file_diag << std::endl;
     std::ofstream out(output_file_diag);
     pcoh.output_diagram(out);
     out.close();
@@ -266,7 +263,7 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool
       "cuboid-file,c", po::value<std::string>(&cuboid_file),
       "Name of file describing the periodic domain. Format is:\n  min_hx min_hy min_hz\n  max_hx max_hy max_hz")(
       "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::cout")(
+      "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)
           ->default_value(std::numeric_limits<Filtration_value>::infinity()),
@@ -288,18 +285,18 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool
   po::notify(vm);
 
   if (vm.count("help") || !vm.count("input-file") || !vm.count("weight-file")) {
-    std::cout << std::endl;
-    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
-    std::cout << "of a 3D Alpha complex defined on a set of input points.\n";
-    std::cout << "3D Alpha complex can be safe (by default) exact or fast, weighted and/or periodic\n\n";
-    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
-    std::cout << "   p   dim b d \n";
-    std::cout << "where dim is the dimension of the homological feature,\n";
-    std::cout << "b and d are respectively the birth and death of the feature and \n";
-    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients.\n\n";
+    std::clog << std::endl;
+    std::clog << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::clog << "of a 3D Alpha complex defined on a set of input points.\n";
+    std::clog << "3D Alpha complex can be safe (by default) exact or fast, weighted and/or periodic\n\n";
+    std::clog << "The output diagram contains one bar per line, written with the convention: \n";
+    std::clog << "   p   dim b d \n";
+    std::clog << "where dim is the dimension of the homological feature,\n";
+    std::clog << "b and d are respectively the birth and death of the feature and \n";
+    std::clog << "p is the characteristic of the field Z/pZ used for homology coefficients.\n\n";
 
-    std::cout << "Usage: " << argv[0] << " [options] input-file weight-file\n\n";
-    std::cout << visible << std::endl;
+    std::clog << "Usage: " << argv[0] << " [options] input-file weight-file\n\n";
+    std::clog << visible << std::endl;
     exit(-1);
   }
 }
diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
index 486347cc..7c898dfd 100644
--- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
@@ -72,13 +72,10 @@ int main(int argc, char **argv) {
   // ----------------------------------------------------------------------------
   // Display information about the alpha complex
   // ----------------------------------------------------------------------------
-  std::cout << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
+  std::clog << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices()
             << " simplices - " << simplex.num_vertices() << " vertices." << std::endl;
 
-  // Sort the simplices in the order of the filtration
-  simplex.initialize_filtration();
-
-  std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl;
+  std::clog << "Simplex_tree dim: " << simplex.dimension() << std::endl;
   // Compute the persistence diagram of the complex
   Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(
       simplex);
@@ -91,7 +88,7 @@ int main(int argc, char **argv) {
   if (output_file_diag.empty()) {
     pcoh.output_diagram();
   } else {
-    std::cout << "Result in file: " << output_file_diag << std::endl;
+    std::clog << "Result in file: " << output_file_diag << std::endl;
     std::ofstream out(output_file_diag);
     pcoh.output_diagram(out);
     out.close();
@@ -114,7 +111,7 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool
       "fast,f", po::bool_switch(&fast),
       "To activate fast version of Alpha complex (default is false, not available if exact is set)")(
       "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::cout")(
+      "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)
                                       ->default_value(std::numeric_limits<Filtration_value>::infinity()),
       "Maximal alpha square value for the Alpha complex construction.")(
@@ -135,17 +132,17 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool
   po::notify(vm);
 
   if (vm.count("help") || !vm.count("input-file")) {
-    std::cout << std::endl;
-    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
-    std::cout << "of an Alpha complex defined on a set of input points.\n \n";
-    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
-    std::cout << "   p   dim b d \n";
-    std::cout << "where dim is the dimension of the homological feature,\n";
-    std::cout << "b and d are respectively the birth and death of the feature and \n";
-    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
-
-    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
-    std::cout << visible << std::endl;
+    std::clog << std::endl;
+    std::clog << "Compute the persistent homology with coefficient field Z/pZ \n";
+    std::clog << "of an Alpha complex defined on a set of input points.\n \n";
+    std::clog << "The output diagram contains one bar per line, written with the convention: \n";
+    std::clog << "   p   dim b d \n";
+    std::clog << "where dim is the dimension of the homological feature,\n";
+    std::clog << "b and d are respectively the birth and death of the feature and \n";
+    std::clog << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::clog << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::clog << visible << std::endl;
     exit(-1);
   }
 }