Use exact option when computing edges filtrations as well in cech

4 files changed, 33 insertions, 14 deletions

diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
index a9dc5d0d..a0e727be 100644
--- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
+++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
@@ -40,9 +40,9 @@ Simplex_tree benchmark_cech(const std::string& off_file_points, const Filtration
     Points_off_reader_cgal off_reader_cgal(off_file_points);
 
     Gudhi::Clock cech_clock("Cech computation");
-    Cech_complex cech_complex_from_points(off_reader_cgal.get_point_cloud(), radius);
+    Cech_complex cech_complex_from_points(off_reader_cgal.get_point_cloud(), radius, exact);
     Simplex_tree cech_stree;
-    cech_complex_from_points.create_complex(cech_stree, dim_max, exact);
+    cech_complex_from_points.create_complex(cech_stree, dim_max);
 
     // ------------------------------------------
     // Display information about the Cech complex
diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h
index 08b7a72f..625f7c9c 100644
--- a/src/Cech_complex/include/gudhi/Cech_complex.h
+++ b/src/Cech_complex/include/gudhi/Cech_complex.h
@@ -62,15 +62,18 @@ class Cech_complex {
    *
    * @param[in] points Range of points where each point is defined as `kernel::Point_d`.
    * @param[in] max_radius Maximal radius value.
+   * @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>.
+   * Default is false.
    *
    */
   template<typename InputPointRange >
-  Cech_complex(const InputPointRange & points, Filtration_value max_radius) : max_radius_(max_radius) {
+  Cech_complex(const InputPointRange & points, Filtration_value max_radius, const bool exact = false) : max_radius_(max_radius), exact_(exact) {
 
     point_cloud_.assign(std::begin(points), std::end(points));
 
     cech_skeleton_graph_ = Gudhi::compute_proximity_graph<SimplicialComplexForCechComplex>(
-        point_cloud_, max_radius_, Sphere_circumradius<Kernel, Filtration_value>());
+        point_cloud_, max_radius_, Sphere_circumradius<Kernel, Filtration_value>(exact));
   }
 
   /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal
@@ -78,19 +81,17 @@ 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, const bool exact = false) {
+  void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) {
     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, exact));
+    complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this));
   }
 
   /** @return max_radius value given at construction. */
@@ -106,11 +107,17 @@ class Cech_complex {
    */
   std::vector<Sphere> & get_cache() { return cache_; }
 
+  /** \brief Check exact option
+   * @return Exact option.
+   */
+  const bool is_exact() { return exact_; }
+
  private:
   Proximity_graph cech_skeleton_graph_;
   Filtration_value max_radius_;
   Point_cloud point_cloud_;
   std::vector<Sphere> cache_;
+  const bool exact_;
 };
 
 }  // namespace cech_complex
diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
index e7f548ba..961da1fd 100644
--- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
+++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
@@ -104,7 +104,7 @@ class Cech_blocker {
 #endif  // DEBUG_TRACES
             is_min_enclos_ball = true;
 #if CGAL_VERSION_NR >= 1050000000
-            if(exact_) CGAL::exact(sph.second);
+            if(cc_ptr_->is_exact()) CGAL::exact(sph.second);
 #endif
             radius = std::sqrt(cast_to_fv(sph.second));
             sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size());
@@ -119,7 +119,7 @@ class Cech_blocker {
         }
         Sphere sph = get_sphere(points.cbegin(), points.cend());
 #if CGAL_VERSION_NR >= 1050000000
-        if(exact_) CGAL::exact(sph.second);
+        if(cc_ptr_->is_exact()) CGAL::exact(sph.second);
 #endif
         radius = std::sqrt(cast_to_fv(sph.second));
 
@@ -135,13 +135,12 @@ class Cech_blocker {
   }
 
   /** \internal \brief Čech complex blocker constructor. */
-  Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr, const bool exact) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr), exact_(exact) {}
+  Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr) {}
 
  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 790f6950..2f916c0a 100644
--- a/src/Cech_complex/include/gudhi/Sphere_circumradius.h
+++ b/src/Cech_complex/include/gudhi/Sphere_circumradius.h
@@ -12,6 +12,7 @@
 #define SPHERE_CIRCUMRADIUS_H_
 
 #include <CGAL/Epick_d.h> // for #include <CGAL/NT_converter.h> which is not working/compiling alone
+#include <CGAL/Lazy_exact_nt.h> // for CGAL::exact
 
 #include <cmath>  // for std::sqrt
 #include <vector>
@@ -26,6 +27,7 @@ template<typename Kernel, typename Filtration_value>
 class Sphere_circumradius {
  private:
     Kernel kernel_;
+    const bool exact_;
  public:
     using FT = typename Kernel::FT;
     using Point = typename Kernel::Point_d;
@@ -42,7 +44,9 @@ class Sphere_circumradius {
    *
    */
   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.;
+    auto squared_dist_obj = kernel_.squared_distance_d_object()(point_1, point_2);
+    if(exact_) CGAL::exact(squared_dist_obj);
+    return std::sqrt(cast_to_fv(squared_dist_obj)) / 2.;
   }
 
   /** \brief Circumradius of sphere passing through point cloud using CGAL.
@@ -53,9 +57,18 @@ class Sphere_circumradius {
    *
    */
   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())));
+    auto squared_radius_obj = kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end());
+    if(exact_) CGAL::exact(squared_radius_obj);
+    return std::sqrt(cast_to_fv(squared_radius_obj));
   }
 
+  /** \brief Constructor
+   * @param[in] exact Option for 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>.
+   * Default is false.
+   */
+  Sphere_circumradius(const bool exact = false) : exact_(exact) {}
+
 };
 
 }  // namespace cech_complex