Merge remote-tracking branch 'upstream/master' into cech_extra_point

5 files changed, 52 insertions, 31 deletions

diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
index d2a71879..a9dc5d0d 100644
--- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
+++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
@@ -31,7 +31,7 @@ using Points_off_reader = Gudhi::Points_off_reader<Point>;
 using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
 
 template<typename Kernel>
-Simplex_tree benchmark_cech(const std::string& off_file_points, const Filtration_value& radius, const int& dim_max) {
+Simplex_tree benchmark_cech(const std::string& off_file_points, const Filtration_value& radius, const int& dim_max, const bool exact) {
     using Point_cgal = typename Kernel::Point_d;
     using Points_off_reader_cgal = Gudhi::Points_off_reader<Point_cgal>;
     using Cech_complex = Gudhi::cech_complex::Cech_complex<Kernel, Simplex_tree>;
@@ -42,7 +42,7 @@ Simplex_tree benchmark_cech(const std::string& off_file_points, const Filtration
     Gudhi::Clock cech_clock("Cech computation");
     Cech_complex cech_complex_from_points(off_reader_cgal.get_point_cloud(), radius);
     Simplex_tree cech_stree;
-    cech_complex_from_points.create_complex(cech_stree, dim_max);
+    cech_complex_from_points.create_complex(cech_stree, dim_max, exact);
 
     // ------------------------------------------
     // Display information about the Cech complex
@@ -56,24 +56,27 @@ int main(int argc, char* argv[]) {
     boost::filesystem::path full_path(boost::filesystem::current_path());
     std::clog << "Current path is : " << full_path << std::endl;
 
-    std::clog << "File name ; Radius ; Rips time ; Dim-3 Epick Cech time ; Dynamic_dim Epick Cech time ; "
-                 "Dim-3 Epeck Cech time ; Dynamic_dim Epeck Cech time ; Cech nb simplices ; Rips nb simplices;"
+    std::clog << "File name ; Radius ; Rips time ; Dim-3 Fast Cech time ; Dynamic_dim Fast Cech time ; "
+                 "Dim-3 Safe Cech time ; Dynamic_dim Safe Cech time ; Dim-3 Exact Cech time ; Dynamic_dim Exact Cech time ; "
+                 "Cech nb simplices ; Rips nb simplices;"
               << std::endl;
     boost::filesystem::directory_iterator end_itr;  // default construction yields past-the-end
+    // For every ".off" file in the current directory, and for 3 predefined thresholds, compare Rips and various Cech constructions
     for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) {
         if (!boost::filesystem::is_directory(itr->status())) {
             if (itr->path().extension() == ".off") {
                 Points_off_reader off_reader(itr->path().string());
                 Point p0 = off_reader.get_point_cloud()[0];
-
-                for (Filtration_value radius = 0.1; radius < 0.4; radius += 0.1) {
+                // Loop over the different thresholds
+                for (Filtration_value radius = 0.1; radius < 0.35; radius += 0.1) {
                     std::clog << itr->path().stem() << "  ;  ";
                     std::clog << radius << "  ;  ";
 
                     Gudhi::Clock rips_clock("Rips computation");
                     Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius, Gudhi::Euclidean_distance());
                     Simplex_tree rips_stree;
-                    rips_complex_from_points.create_complex(rips_stree, p0.size() - 1);
+                    int dim_max = p0.size() - 1;
+                    rips_complex_from_points.create_complex(rips_stree, dim_max);
                     // ------------------------------------------
                     // Display information about the Rips complex
                     // ------------------------------------------
@@ -83,10 +86,15 @@ int main(int argc, char* argv[]) {
                     // --------------
                     // Cech complex
                     // --------------
-                    benchmark_cech<CGAL::Epick_d<CGAL::Dimension_tag<3>>>(itr->path().string(), radius, p0.size() - 1);
-                    benchmark_cech<CGAL::Epick_d<CGAL::Dynamic_dimension_tag>>(itr->path().string(), radius, p0.size() - 1);
-                    benchmark_cech<CGAL::Epeck_d<CGAL::Dimension_tag<3>>>(itr->path().string(), radius, p0.size() - 1);
-                    auto cech_stree = benchmark_cech<CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>(itr->path().string(), radius, p0.size() - 1);
+                    // Fast
+                    benchmark_cech<CGAL::Epick_d<CGAL::Dimension_tag<3>>>(itr->path().string(), radius, dim_max, false);
+                    benchmark_cech<CGAL::Epick_d<CGAL::Dynamic_dimension_tag>>(itr->path().string(), radius, dim_max, false);
+                    // Safe
+                    benchmark_cech<CGAL::Epeck_d<CGAL::Dimension_tag<3>>>(itr->path().string(), radius, dim_max, false);
+                    benchmark_cech<CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>(itr->path().string(), radius, dim_max, false);
+                    // Exact
+                    benchmark_cech<CGAL::Epeck_d<CGAL::Dimension_tag<3>>>(itr->path().string(), radius, dim_max, true);
+                    auto cech_stree = benchmark_cech<CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>(itr->path().string(), radius, dim_max, true);
 
                     std::clog << cech_stree.num_simplices() << "  ;  ";
                     std::clog << rips_stree.num_simplices() << ";" << std::endl;
diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h
index fc39f75b..08b7a72f 100644
--- a/src/Cech_complex/include/gudhi/Cech_complex.h
+++ b/src/Cech_complex/include/gudhi/Cech_complex.h
@@ -30,7 +30,7 @@ namespace cech_complex {
  * \ingroup cech_complex
  *
  * \details
- * Cech complex is a simplicial complex constructed from a proximity graph, where the set of all simplices is filtered
+ * Cech complex is a simplicial complex where the set of all simplices is filtered
  * by the radius of their minimal enclosing ball and bounded by the given max_radius.
  *
  * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d.
@@ -70,7 +70,7 @@ class Cech_complex {
     point_cloud_.assign(std::begin(points), std::end(points));
 
     cech_skeleton_graph_ = Gudhi::compute_proximity_graph<SimplicialComplexForCechComplex>(
-        point_cloud_, max_radius_, Sphere_circumradius<Kernel>());
+        point_cloud_, max_radius_, Sphere_circumradius<Kernel, Filtration_value>());
   }
 
   /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal
@@ -78,17 +78,19 @@ class Cech_complex {
    *
    * @param[in] complex SimplicialComplexForCech to be created.
    * @param[in] dim_max graph expansion until this given maximal dimension.
+   * @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>.
    * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0.
    *
    */
-  void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) {
+  void create_complex(SimplicialComplexForCechComplex& complex, int dim_max, const bool exact = false) {
     GUDHI_CHECK(complex.num_vertices() == 0,
                 std::invalid_argument("Cech_complex::create_complex - simplicial complex is not empty"));
 
     // insert the proximity graph in the simplicial complex
     complex.insert_graph(cech_skeleton_graph_);
     // expand the graph until dimension dim_max
-    complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this));
+    complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this, exact));
   }
 
   /** @return max_radius value given at construction. */
diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
index 9917999f..1bb205b3 100644
--- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
+++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
@@ -12,6 +12,7 @@
 #define CECH_COMPLEX_BLOCKER_H_
 
 #include <CGAL/NT_converter.h> // for casting from FT to Filtration_value
+#include <CGAL/Lazy_exact_nt.h> // for CGAL::exact
 
 #include <iostream>
 #include <vector>
@@ -84,9 +85,9 @@ class Cech_blocker {
             Point_cloud face_points;
             for (auto vertex : sc_ptr_->simplex_vertex_range(face_opposite_vertex.first)) {
                 face_points.push_back(cc_ptr_->get_point(vertex));
-                #ifdef DEBUG_TRACES
-                    std::clog << "#(" << vertex << ")#";
-                #endif  // DEBUG_TRACES
+#ifdef DEBUG_TRACES
+                std::clog << "#(" << vertex << ")#";
+#endif  // DEBUG_TRACES
             }
             sph = get_sphere(face_points.cbegin(), face_points.cend());
             // Put edge sphere in cache
@@ -97,10 +98,13 @@ class Cech_blocker {
         }
         // Check if the minimal enclosing ball of current face contains the extra point/opposite vertex
         if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(face_opposite_vertex.second)) <= sph.second) {
-            #ifdef DEBUG_TRACES
-                std::clog << "center: " << sph.first << ", radius: " <<  radius << std::endl;
-            #endif  // DEBUG_TRACES
+#ifdef DEBUG_TRACES
+            std::clog << "center: " << sph.first << ", radius: " <<  radius << std::endl;
+#endif  // DEBUG_TRACES
             is_min_enclos_ball = true;
+#if CGAL_VERSION_NR >= 1050000000
+            if(exact_) CGAL::exact(sph.second);
+#endif
             radius = std::sqrt(cast_to_fv(sph.second));
             sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size());
             cc_ptr_->get_cache().push_back(sph);
@@ -114,10 +118,13 @@ class Cech_blocker {
             points.push_back(cc_ptr_->get_point(vertex));
         }
         Sphere sph = get_sphere(points.cbegin(), points.cend());
+#if CGAL_VERSION_NR >= 1050000000
+        if(exact_) CGAL::exact(sph.second);
+#endif
         radius = std::sqrt(cast_to_fv(sph.second));
 
         sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size());
-        cc_ptr_->get_cache().push_back(sph);
+        cc_ptr_->get_cache().push_back(std::move(sph));
     }
 
 #ifdef DEBUG_TRACES
@@ -128,12 +135,13 @@ class Cech_blocker {
   }
 
   /** \internal \brief Čech complex blocker constructor. */
-  Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr) {}
+  Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr, const bool exact) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr), exact_(exact) {}
 
  private:
   SimplicialComplexForCech* sc_ptr_;
   Cech_complex* cc_ptr_;
   Kernel kernel_;
+  const bool exact_;
 };
 
 }  // namespace cech_complex
diff --git a/src/Cech_complex/include/gudhi/Sphere_circumradius.h b/src/Cech_complex/include/gudhi/Sphere_circumradius.h
index b0d9f7cc..790f6950 100644
--- a/src/Cech_complex/include/gudhi/Sphere_circumradius.h
+++ b/src/Cech_complex/include/gudhi/Sphere_circumradius.h
@@ -11,7 +11,7 @@
 #ifndef SPHERE_CIRCUMRADIUS_H_
 #define SPHERE_CIRCUMRADIUS_H_
 
-#include <CGAL/Epeck_d.h> // for #include <CGAL/NewKernel_d/KernelD_converter.h>
+#include <CGAL/Epick_d.h> // for #include <CGAL/NT_converter.h> which is not working/compiling alone
 
 #include <cmath>  // for std::sqrt
 #include <vector>
@@ -22,14 +22,17 @@ namespace cech_complex {
 
 /** \private @brief Compute the circumradius of the sphere passing through points given by a range of coordinates.
  * The points are assumed to have the same dimension. */
-template<typename Kernel>
+template<typename Kernel, typename Filtration_value>
 class Sphere_circumradius {
  private:
     Kernel kernel_;
  public:
+    using FT = typename Kernel::FT;
     using Point = typename Kernel::Point_d;
     using Point_cloud = typename std::vector<Point>;
 
+    CGAL::NT_converter<FT, Filtration_value> cast_to_fv;
+
    /** \brief Circumradius of sphere passing through two points using CGAL.
    *
    * @param[in] point_1
@@ -38,8 +41,8 @@ class Sphere_circumradius {
    * \tparam Point must be a Kernel::Point_d from CGAL.
    *
    */
-  double operator()(const Point& point_1, const Point& point_2) const {
-    return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.;
+  Filtration_value operator()(const Point& point_1, const Point& point_2) const {
+    return std::sqrt(cast_to_fv(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.;
   }
 
   /** \brief Circumradius of sphere passing through point cloud using CGAL.
@@ -49,8 +52,8 @@ class Sphere_circumradius {
    * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL.
    *
    */
-  double operator()(const Point_cloud& point_cloud) const {
-    return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end())));
+  Filtration_value operator()(const Point_cloud& point_cloud) const {
+    return std::sqrt(cast_to_fv(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end())));
   }
 
 };
diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp
index ea32f596..f5980e6d 100644
--- a/src/Cech_complex/test/test_cech_complex.cpp
+++ b/src/Cech_complex/test/test_cech_complex.cpp
@@ -107,11 +107,11 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) {
         std::clog << vertex << ",";
         vp.push_back(points.at(vertex));
       }
-      std::clog << ") - distance =" << Gudhi::cech_complex::Sphere_circumradius<Kernel>()(vp.at(0), vp.at(1))
+      std::clog << ") - distance =" << Gudhi::cech_complex::Sphere_circumradius<Kernel, Filtration_value>()(vp.at(0), vp.at(1))
                 << " - filtration =" << st.filtration(f_simplex) << std::endl;
       BOOST_CHECK(vp.size() == 2);
       GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex),
-                                      Gudhi::cech_complex::Sphere_circumradius<Kernel>()(vp.at(0), vp.at(1)));
+                                      Gudhi::cech_complex::Sphere_circumradius<Kernel, Filtration_value>()(vp.at(0), vp.at(1)));
     }
   }