Modify cech benchmark to include both Epick and Epeck cases

1 files changed, 109 insertions, 68 deletions

diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
index bf013a81..9cf24542 100644
--- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
+++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
@@ -18,6 +18,7 @@
 #include <gudhi/Miniball.hpp>
 
 #include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h>
 
 #include "boost/filesystem.hpp"  // includes all needed Boost.Filesystem declarations
 
@@ -32,10 +33,6 @@ using Point_cloud = std::vector<Point>;
 using Points_off_reader = Gudhi::Points_off_reader<Point>;
 using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>;
 using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
-using Kernel =  CGAL::Epick_d<CGAL::Dimension_tag<3>>;
-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>;
 
 class Minimal_enclosing_ball_radius {
  public:
@@ -61,79 +58,123 @@ class Minimal_enclosing_ball_radius {
   }
 };
 
-int main(int argc, char* argv[]) {
-  std::string off_file_points = "tore3D_1307.off";
-  Filtration_value threshold = 1e20;
+enum distance_type { Euclidean_dist, Minimal_enclosing_ball_dist, CGAL_dist };
 
-  // Extract the points from the file filepoints
-  Points_off_reader off_reader(off_file_points);
-  Points_off_reader_cgal off_reader_cgal(off_file_points);
+template<bool use_cgal, typename Kernel = std::enable_if_t<use_cgal>>
+void benchmark_prox_graph(const std::string& off_file_points, const Filtration_value& threshold, const std::string& msg, distance_type dist = CGAL_dist) {
+    if (dist != CGAL_dist) {
+        std::cerr << "Error: when CGAL is used, the distance should be CGAL_dist" << std::endl;
+        exit(-1);
+    }
+    if (!use_cgal) {
+        std::cerr << "Warning: if kernel is given, CGAL will be used" << std::endl;
+    }
+    using Point_cgal = typename Kernel::Point_d;
+    using Points_off_reader_cgal = Gudhi::Points_off_reader<Point_cgal>;
 
-  Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance");
-  // Compute the proximity graph of the points
-  Proximity_graph euclidean_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(
-      off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance());
+    // Extract the points from the file filepoints
+    Points_off_reader_cgal off_reader_cgal(off_file_points);
 
-  std::clog << euclidean_clock << std::endl;
+    Gudhi::Clock cgal_circumsphere_clock("Gudhi::cech_complex::Sphere_circumradius_cgal()");
+    // Compute the proximity graph of the points
+    Proximity_graph cgal_circumsphere_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader_cgal.get_point_cloud(), threshold,
+                                                                                                Gudhi::cech_complex::Sphere_circumradius<Kernel>());
+    std::clog << msg << " - " << cgal_circumsphere_clock << std::endl;
+}
 
-  Gudhi::Clock miniball_clock("Minimal_enclosing_ball_radius");
-  // Compute the proximity graph of the points
-  Proximity_graph miniball_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(
-      off_reader.get_point_cloud(), threshold, Minimal_enclosing_ball_radius());
-  std::clog << miniball_clock << std::endl;
+template<bool use_cgal>
+void benchmark_prox_graph(const std::string& off_file_points, const Filtration_value& threshold, const std::string& msg, distance_type dist) {
+    // Extract the points from the file filepoints
+    Points_off_reader off_reader(off_file_points);
+
+    if (dist == Euclidean_dist) {
+        Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance");
+        // Compute the proximity graph of the points
+        Proximity_graph euclidean_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), threshold,
+                                                                                                    Gudhi::Euclidean_distance());
+        std::clog << msg << " - " << euclidean_clock << std::endl;
+    }
+    else if (dist == Minimal_enclosing_ball_dist) {
+        Gudhi::Clock miniball_clock("Minimal_enclosing_ball_radius");
+        // Compute the proximity graph of the points
+        Proximity_graph miniball_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), threshold,
+                                                                                           Minimal_enclosing_ball_radius());
+        std::clog << msg << " - " << miniball_clock << std::endl;
+    }
+    else {
+        std::cerr << "Error: when CGAL is not used, the distance should be either Euclidean_dist or Minimal_enclosing_ball_dist" << std::endl;
+        exit(-1);
+    }
+}
+
+template<typename Kernel>
+void benchmark_cech(const std::string& off_file_points, const Filtration_value& radius, const int& dim_max) {
+    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>;
+
+    // Extract the points from the file filepoints
+    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);
+    Simplex_tree cech_stree;
+    cech_complex_from_points.create_complex(cech_stree, dim_max);
+
+    // ------------------------------------------
+    // Display information about the Cech complex
+    // ------------------------------------------
+    double cech_sec = cech_clock.num_seconds();
+    std::clog << cech_sec << "  ;  ";
+    std::clog << cech_stree.num_simplices() << "  ;  ";
+}
 
-  Gudhi::Clock cgal_circumsphere_clock("Gudhi::cech_complex::Sphere_circumradius_cgal()");
-  // Compute the proximity graph of the points
-  Proximity_graph cgal_circumsphere_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(
-      off_reader_cgal.get_point_cloud(), threshold, Gudhi::cech_complex::Sphere_circumradius<Kernel>());
-  std::clog << cgal_circumsphere_clock << std::endl;
+int main(int argc, char* argv[]) {
+    std::string off_file_points = "tore3D_1307.off";
+    Filtration_value threshold = 1e20;
+
+    benchmark_prox_graph<false>(off_file_points, threshold, "Euclidean distance", Euclidean_dist);
+    benchmark_prox_graph<false>(off_file_points, threshold, "Minimal_enclosing_ball", Minimal_enclosing_ball_dist);
+    benchmark_prox_graph<true, CGAL::Epick_d<CGAL::Dimension_tag<3>>>(off_file_points, threshold, "Epick");
+    benchmark_prox_graph<true, CGAL::Epeck_d<CGAL::Dimension_tag<3>>>(off_file_points, threshold, "Epeck");
 
-  boost::filesystem::path full_path(boost::filesystem::current_path());
-  std::clog << "Current path is : " << full_path << std::endl;
+    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; Cech time; Ratio Rips/Cech time; Rips nb simplices; Cech nb simplices;"
+    std::clog << "File name ; Radius ; Rips time ; Epick Cech time ; Epick Cech nb simplices ; Epeck Cech time ; Epeck Cech nb simplices ; Rips nb simplices;"
             << std::endl;
-  boost::filesystem::directory_iterator end_itr;  // default construction yields past-the-end
-  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")  // see below
-      {
-        Points_off_reader off_reader(itr->path().string());
-        Points_off_reader_cgal off_reader_cgal(itr->path().string());
-
-        Point p0 = off_reader.get_point_cloud()[0];
-
-        for (Filtration_value radius = 0.1; radius < 0.4; 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);
-          // ------------------------------------------
-          // Display information about the Rips complex
-          // ------------------------------------------
-          double rips_sec = rips_clock.num_seconds();
-          std::clog << rips_sec << ";";
-
-          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, p0.size() - 1);
-          // ------------------------------------------
-          // Display information about the Cech complex
-          // ------------------------------------------
-          double cech_sec = cech_clock.num_seconds();
-          std::clog << cech_sec << ";";
-          std::clog << cech_sec / rips_sec << ";";
-
-          assert(rips_stree.num_simplices() >= cech_stree.num_simplices());
-          std::clog << rips_stree.num_simplices() << ";";
-          std::clog << cech_stree.num_simplices() << ";" << std::endl;
+    boost::filesystem::directory_iterator end_itr;  // default construction yields past-the-end
+    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) {
+                    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);
+                    // ------------------------------------------
+                    // Display information about the Rips complex
+                    // ------------------------------------------
+                    double rips_sec = rips_clock.num_seconds();
+                    std::clog << rips_sec << "  ;  ";
+
+                    // --------------
+                    // Cech complex
+                    // --------------
+                    benchmark_cech<CGAL::Epick_d<CGAL::Dimension_tag<3>>>(itr->path().string(), radius, p0.size() - 1);
+                    benchmark_cech<CGAL::Epeck_d<CGAL::Dimension_tag<3>>>(itr->path().string(), radius, p0.size() - 1);
+
+                    std::clog << rips_stree.num_simplices() << ";" << std::endl;
+                }
+            }
         }
-      }
     }
-  }
 
-  return 0;
+    return 0;
 }