Move minimal enclosing balls cache to Cech_complex.h instead of the blocker

Modify cech test to be more relevant regarding the current algorithm
Do some cleaning

10 files changed, 54 insertions, 225 deletions

diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
index c332c656..4a1aa06e 100644
--- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
+++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
@@ -33,7 +33,7 @@ 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::Epeck_d<CGAL::Dynamic_dimension_tag>;
-using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel>;
+using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel, Simplex_tree>;
 
 class Minimal_enclosing_ball_radius {
  public:
diff --git a/src/Cech_complex/example/cech_complex_example_from_points.cpp b/src/Cech_complex/example/cech_complex_example_from_points.cpp
index e78ad51d..78861951 100644
--- a/src/Cech_complex/example/cech_complex_example_from_points.cpp
+++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp
@@ -10,25 +10,15 @@
 
 int main() {
   // Type definitions
-//   using Point_cloud = std::vector<std::array<double, 2>>;
   using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
   using Filtration_value = Simplex_tree::Filtration_value;
   using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
   using FT = typename Kernel::FT;
   using Point = typename Kernel::Point_d;
   using Point_cloud = std::vector<Point>;
-  using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel>;
+  using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel, Simplex_tree>;
 
   Point_cloud points;
-//   points.push_back({1., 0.});                  // 0
-//   points.push_back({0., 1.});                  // 1
-//   points.push_back({2., 1.});                  // 2
-//   points.push_back({3., 2.});                  // 3
-//   points.push_back({0., 3.});                  // 4
-//   points.push_back({3. + std::sqrt(3.), 3.});  // 5
-
-//   std::vector<FT> point({0.0, 0.0, 0.0, 0.0});
-//   points.emplace_back(point.begin(), point.end());
 
   std::vector<FT> point0({1., 0.});
   points.emplace_back(point0.begin(), point0.end());
@@ -42,8 +32,6 @@ int main() {
   points.emplace_back(point4.begin(), point4.end());
   std::vector<FT> point5({3. + std::sqrt(3.), 3.});
   points.emplace_back(point5.begin(), point5.end());
-
-/*
   std::vector<FT> point6({1., 4.});
   points.emplace_back(point6.begin(), point6.end());
   std::vector<FT> point7({3., 4.});
@@ -54,39 +42,15 @@ int main() {
   points.emplace_back(point9.begin(), point9.end());
   std::vector<FT> point10({-0.5, 2.});
   points.emplace_back(point10.begin(), point10.end());
-*/
-//   points.emplace_back(Point(std::vector<FT>({1., 0.})));
-//   points.emplace_back(Point(std::vector<FT>({0., 1.})));
-//   points.emplace_back(Point(std::vector<FT>({2., 1.})));
-//   points.emplace_back(Point(std::vector<FT>({3., 2.})));
-//   points.emplace_back(Point(std::vector<FT>({0., 3.})));
-//   points.emplace_back(Point(std::vector<FT>({3. + std::sqrt(3.), 3.})));
-
-
-//   points.push_back(Point(1.0, 0.0));
-//   points.push_back(Point(0.0, 1.0));
-//   points.push_back(Point(2.0, 1.0));
-//   points.push_back(Point(3.0, 2.0));
-//   points.push_back(Point(0.0, 3.0));
-//   points.push_back(Point(3.0 + std::sqrt(3.0), 3.0));
-
-
-//   points.push_back({1., 4.});                  // 6
-//   points.push_back({3., 4.});                  // 7
-//   points.push_back({2., 4. + std::sqrt(3.)});  // 8
-//   points.push_back({0., 4.});                  // 9
-//   points.push_back({-0.5, 2.});                // 10
 
   // ----------------------------------------------------------------------------
   // Init of a Cech complex from points
   // ----------------------------------------------------------------------------
-  Filtration_value max_radius = 10.; //100.;
-  std::clog << "Hind: Just before the Cech constructor" << std::endl;
+  Filtration_value max_radius = 100.; //100.;
   Cech_complex cech_complex_from_points(points, max_radius);
-  std::clog << "Hind: Just after the Cech constructor" << std::endl;
 
   Simplex_tree stree;
-  cech_complex_from_points.create_complex(stree, 3); //6
+  cech_complex_from_points.create_complex(stree, 6); //6
   // ----------------------------------------------------------------------------
   // Display information about the one skeleton Cech complex
   // ----------------------------------------------------------------------------
diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp
index ac08e6cc..44e7f945 100644
--- a/src/Cech_complex/example/cech_complex_step_by_step.cpp
+++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp
@@ -13,8 +13,6 @@
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Points_off_io.h>
 
-#include <gudhi/Miniball.hpp> // TODO to remove ?
-
 #include <CGAL/Epeck_d.h>
 
 #include <boost/program_options.hpp>
@@ -36,7 +34,6 @@
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Simplex_handle = Simplex_tree::Simplex_handle;
 using Filtration_value = Simplex_tree::Filtration_value;
-// using Point = std::vector<double>;
 using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
 using Point = typename Kernel::Point_d;
 using Points_off_reader = Gudhi::Points_off_reader<Point>;
@@ -45,9 +42,6 @@ using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>;
 class Cech_blocker {
  private:
   using Point_cloud = std::vector<Point>;
-//   using Point_iterator = Point_cloud::const_iterator;
-//   using Coordinate_iterator = Point::const_iterator;
-//   using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
 
  public:
   bool operator()(Simplex_handle sh) {
@@ -67,14 +61,12 @@ class Cech_blocker {
   }
   Cech_blocker(Simplex_tree& simplex_tree, Filtration_value max_radius, const std::vector<Point>& point_cloud)
       : simplex_tree_(simplex_tree), max_radius_(max_radius), point_cloud_(point_cloud) {
-//     dimension_ = point_cloud_[0].size();
   }
 
  private:
   Simplex_tree simplex_tree_;
   Filtration_value max_radius_;
   std::vector<Point> point_cloud_;
-//   int dimension_;
 };
 
 void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& max_radius, int& dim_max);
diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h
index 2c6f202a..32a78aec 100644
--- a/src/Cech_complex/include/gudhi/Cech_complex.h
+++ b/src/Cech_complex/include/gudhi/Cech_complex.h
@@ -40,7 +40,7 @@ namespace cech_complex {
  * \tparam ForwardPointRange must be a range for which `std::begin()` and `std::end()` methods return input
  * iterators on a point. `std::begin()` and `std::end()` methods are also required for a point.
  */
-template <typename SimplicialComplexForProximityGraph, typename ForwardPointRange, typename Kernel>
+template <typename SimplicialComplexForProximityGraph, typename ForwardPointRange, typename Kernel, typename SimplicialComplexForCechComplex>
 class Cech_complex {
  private:
   // Required by compute_proximity_graph
@@ -48,17 +48,17 @@ class Cech_complex {
   using Filtration_value = typename SimplicialComplexForProximityGraph::Filtration_value;
   using Proximity_graph = Gudhi::Proximity_graph<SimplicialComplexForProximityGraph>;
 
-  // Retrieve Coordinate type from ForwardPointRange
-//   using Point_from_range_iterator = typename boost::range_const_iterator<ForwardPointRange>::type;
-//   using Point_from_range = typename std::iterator_traits<Point_from_range_iterator>::value_type;
-//   using Coordinate_iterator = typename boost::range_const_iterator<Point_from_range>::type;
-//   using Coordinate = typename std::iterator_traits<Coordinate_iterator>::value_type;
-
  public:
-  // Point and Point_cloud type definition
-  //using Point = std::vector<Coordinate>;
-  using Point = typename Kernel::Point_d;
-  using Point_cloud = std::vector<Point>;
+
+  using cech_blocker = Cech_blocker<SimplicialComplexForCechComplex, Cech_complex, Kernel>;
+
+  using Point_d = typename cech_blocker::Point_d;
+  using Point_cloud = std::vector<Point_d>;
+
+  // Numeric type of coordinates in the kernel
+  using FT = typename cech_blocker::FT;
+  // Sphere is a pair of point and squared radius.
+  using Sphere = typename std::pair<Point_d, FT>;
 
  public:
   /** \brief Cech_complex constructor from a list of points.
@@ -77,7 +77,6 @@ class Cech_complex {
 
     point_cloud_.assign(points.begin(), points.end());
 
-    std::clog << "Hind: Just before the graph compute" << std::endl;
     cech_skeleton_graph_ = Gudhi::compute_proximity_graph<SimplicialComplexForProximityGraph>(
         point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius());
   }
@@ -90,19 +89,14 @@ class Cech_complex {
    * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0.
    *
    */
-  template <typename SimplicialComplexForCechComplex>
   void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) {
-    std::clog << "Hind: in create complex" << std::endl;
     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
-    std::clog << "Hind: before insert_graph" << std::endl;
     complex.insert_graph(cech_skeleton_graph_);
     // expand the graph until dimension dim_max
-    std::clog << "Hind: before expansion_with_blockers" << std::endl;
-    complex.expansion_with_blockers(dim_max,
-                                    Cech_blocker<SimplicialComplexForCechComplex, Cech_complex, Kernel>(&complex, this));
+    complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this));
   }
 
   /** @return max_radius value given at construction. */
@@ -111,12 +105,18 @@ class Cech_complex {
   /** @param[in] vertex Point position in the range.
    * @return The point.
    */
-  const Point& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; }
+  const Point_d& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; }
+
+  /**
+   * @return Vector of cached spheres.
+   */
+  std::vector<Sphere> & get_cache() { return cache_; }
 
  private:
   Proximity_graph cech_skeleton_graph_;
   Filtration_value max_radius_;
   Point_cloud point_cloud_;
+  std::vector<Sphere> cache_;
 };
 
 }  // 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 3cac9ee2..5edd005d 100644
--- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
+++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
@@ -11,20 +11,11 @@
 #ifndef CECH_COMPLEX_BLOCKER_H_
 #define CECH_COMPLEX_BLOCKER_H_
 
-// TODO to remove
-#include <gudhi/distance_functions.h>  // for Gudhi::Minimal_enclosing_ball_radius
-
-#include <CGAL/number_utils.h> // for CGAL::to_double
-#include <CGAL/NT_converter.h> //
-// #include <CGAL/Delaunay_triangulation.h>
-#include <CGAL/Epeck_d.h>  // For EXACT or SAFE version
-#include <CGAL/Spatial_sort_traits_adapter_d.h>
+#include <CGAL/NT_converter.h> // for casting from FT to double
 
 #include <iostream>
 #include <vector>
 #include <set>
-#include <map>
-#include <algorithm>
 #include <cmath>  // for std::sqrt
 
 namespace Gudhi {
@@ -48,7 +39,6 @@ namespace cech_complex {
 template <typename SimplicialComplexForCech, typename Cech_complex, typename Kernel>
 class Cech_blocker {
  private:
-//   using Point_cloud = typename Cech_complex::Point_cloud;
 
   using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle;
   using Filtration_value = typename SimplicialComplexForCech::Filtration_value;
@@ -61,28 +51,18 @@ class Cech_blocker {
   // Sphere is a pair of point and squared radius.
   using Sphere = typename std::pair<Point_d, FT>;
 
-
-   /** \internal \brief TODO
-   *  \param[in] 
-   *  \return  */
   template<class PointIterator>
   FT get_squared_radius(PointIterator begin, PointIterator end) const {
     return kernel_.compute_squared_radius_d_object()(begin, end);
   }
-  
-  /** \internal \brief TODO
-   *  \param[in] 
-   *  \return  */
+
   template<class PointIterator>
   Sphere get_sphere(PointIterator begin, PointIterator end) const {
     Point_d c = kernel_.construct_circumcenter_d_object()(begin, end);
     FT r = kernel_.squared_distance_d_object()(c, *begin);
     return std::make_pair(std::move(c), std::move(r));
   }
-  
-  /** \internal \brief TODO
-   *  \param[in] 
-   *  \return  */
+
   template<typename Sphere>
   class CompareSpheresRadii
   {
@@ -93,7 +73,7 @@ class Cech_blocker {
         return cast_to_double(firstSphere.second) < cast_to_double(secondSphere.second);
       }
   };
-  
+
   /** \internal \brief Čech complex blocker operator() - the oracle - assigns the filtration value from the simplex
    * radius and returns if the simplex expansion must be blocked.
    *  \param[in] sh The Simplex_handle.
@@ -102,104 +82,54 @@ class Cech_blocker {
     using Point_cloud =  std::vector<Point_d>;
     CGAL::NT_converter<FT, double> cast_to_double;
     Filtration_value radius = 0.;
-//     std::string key_to_permute;
-//     std::vector<std::string> faces_keys;
-    
+
     // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices
-    // std::set <Filtration_value> enclosing_ball_radii;
     std::set <Sphere, CompareSpheresRadii<Sphere>> enclosing_ball_spheres;
     for (auto face : sc_ptr_->boundary_simplex_range(sh)) {
-        
         // Find which vertex of sh is missing in face. We rely on the fact that simplex_vertex_range is sorted.
         auto longlist = sc_ptr_->simplex_vertex_range(sh);
         auto shortlist = sc_ptr_->simplex_vertex_range(face);
 
-//         std::clog << "Hind debug: within FACE loop "<< std::endl;
-        // TODO to remove
-//         for (auto i = std::begin(longlist); i != std::end(longlist);++i)
-//             std::clog << "Hind debug: longlist: " << cc_ptr_->get_point(*i) << std::endl;
-//         for (auto i = std::begin(shortlist); i != std::end(shortlist);++i)
-//             std::clog << "Hind debug: shortlist: " << cc_ptr_->get_point(*i) << std::endl;
-
         auto longiter = std::begin(longlist);
         auto shortiter = std::begin(shortlist);
         auto enditer = std::end(shortlist);
         while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; }
         auto extra = *longiter; // Vertex_handle
 
-//         std::clog << "Hind debug: extra vertex: " << cc_ptr_->get_point(extra) << std::endl;
-        
         Point_cloud face_points;
-//         std::string key, key_extra;
         for (auto vertex : sc_ptr_->simplex_vertex_range(face)) {
             face_points.push_back(cc_ptr_->get_point(vertex));
-//             key.append(std::to_string(vertex));
     #ifdef DEBUG_TRACES
         std::clog << "#(" << vertex << ")#";
     #endif  // DEBUG_TRACES
         }
-//         key_extra = key;
-//         key_extra.append(std::to_string(extra));
-//         faces_keys.push_back(key_extra);
-//         key_to_permute = key_extra;
-//         std::clog << "END OF VERTICES " << std::endl;
-//         std::clog << "KEY is: " << key << std::endl;
-//         std::clog << "KEY extra is: " << key_extra << std::endl;
         Sphere sph;
-        auto k = sc_ptr_->key(sh);
-        if(k != sc_ptr_->null_key())
-            sph = cache_[k];
+        auto face_sh = sc_ptr_->find(sc_ptr_->simplex_vertex_range(face));
+        auto k = sc_ptr_->key(face_sh);
+        if(k != sc_ptr_->null_key()) {
+            sph = cc_ptr_->get_cache().at(k);
+        }
         else {
             sph = get_sphere(face_points.cbegin(), face_points.cend());
         }
-//         auto it = cache_.find(key);
-//         if(it != cache_.end())
-//             sph = it->second;
-//         else {
-//             sph = get_sphere(face_points.cbegin(), face_points.cend());
-//         }
 
         if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) {
             radius = std::sqrt(cast_to_double(sph.second));
             #ifdef DEBUG_TRACES
                 std::clog << "circumcenter: " << sph.first << ", radius: " <<  radius << std::endl;
             #endif  // DEBUG_TRACES
-//             std::clog << "distance FYI: " << kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " < " << cast_to_double(sph.second) << std::endl;
-            // enclosing_ball_radii.insert(radius);
             enclosing_ball_spheres.insert(sph);
-//             cache_[key_extra] = sph;
-//             sc_ptr_->assign_key(sh, cache_.size());
-//             cache_.push_back(sph);
         }
-//         else {// TODO to remove
-//             std::clog << "vertex not included BECAUSE DISTANCE: "<< kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " AND RAD SPHERE: " << sph.second << std::endl;
-//         }
     }
     // Get the minimal radius of all faces enclosing balls if exists
     if (!enclosing_ball_spheres.empty()) {
-        // radius = *enclosing_ball_radii.begin();
         Sphere sph_min = *enclosing_ball_spheres.begin();
         radius = std::sqrt(cast_to_double(sph_min.second));
-        // std::clog << "CHECK that radius of min sphere is min radius: " << std::sqrt(cast_to_double(sph_min.second)) << "; and RADIUS min: " << radius << std::endl;
-        // Set all key_to_permute permutations to min sphere in cache
-//         do
-//         {
-//             if (cache_.find(key_to_permute) != cache_.end()) {
-//                 if (cast_to_double(cache_[key_to_permute].second) >  cast_to_double(sph_min.second))
-//                     cache_[key_to_permute] = sph_min;
-//             }
-//             else {
-//                 cache_[key_to_permute] = sph_min;
-//             }
-//         } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end()));
-//         for (auto k : faces_keys) {
-//             cache_[k] = sph_min;
-//         }
-        sc_ptr_->assign_key(sh, cache_.size());
-        cache_.push_back(sph_min);
+
+        sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size());
+        cc_ptr_->get_cache().push_back(sph_min);
     }
-//     std::clog << "END OF FACES ; radius = " << radius << std::endl;
-    
+
     if (radius == 0.) { // Spheres of each face don't contain the whole simplex
         Point_cloud points;
         for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) {
@@ -207,25 +137,10 @@ class Cech_blocker {
         }
         Sphere sph = get_sphere(points.cbegin(), points.cend());
         radius = std::sqrt(cast_to_double(sph.second));
-//         std::clog << "GLOBAL SPHERE radius = " << radius << std::endl;
-        // Set all key_to_permute permutations to sphere in cache
-//         do
-//         {
-// //             if (cache_.find(key_to_permute) != cache_.end()) {
-// //                 if (cast_to_double(cache_[key_to_permute].second) >  cast_to_double(sph.second))
-// //                     cache_[key_to_permute] = sph;
-// //             }
-// //             else {
-// //                 cache_[key_to_permute] = sph;
-// //             }
-//         } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end()));
-//         for (auto k : faces_keys) {
-//             cache_[k] = sph;
-//         }
-//         sc_ptr_->assign_key(sh, cache_.size());
-//         cache_.push_back(sph);
-    }
 
+        sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size());
+        cc_ptr_->get_cache().push_back(sph);
+    }
 
 #ifdef DEBUG_TRACES
     if (radius > cc_ptr_->max_radius()) std::clog << "radius > max_radius => expansion is blocked\n";
@@ -241,8 +156,6 @@ class Cech_blocker {
   SimplicialComplexForCech* sc_ptr_;
   Cech_complex* cc_ptr_;
   Kernel kernel_;
-//   std::map<std::string, Sphere> cache_;
-  std::vector<Sphere> cache_;
 };
 
 }  // namespace cech_complex
diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp
index 81efd6ae..51b466da 100644
--- a/src/Cech_complex/test/test_cech_complex.cpp
+++ b/src/Cech_complex/test/test_cech_complex.cpp
@@ -24,26 +24,19 @@
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/distance_functions.h>
 #include <gudhi/Unitary_tests_utils.h>
-#include <gudhi/Miniball.hpp>
 
 #include <CGAL/Epeck_d.h>  // For EXACT or SAFE version
 
-
 // Type definitions
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Filtration_value = Simplex_tree::Filtration_value;
-//using Point = std::vector<Filtration_value>;
 using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
 using FT = typename Kernel::FT;
 using Point = typename Kernel::Point_d;
 
 using Point_cloud = std::vector<Point>;
 using Points_off_reader = Gudhi::Points_off_reader<Point>;
-using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel>;
-
-// using Point_iterator = Point_cloud::const_iterator;
-// using Coordinate_iterator = Point::const_iterator;
-// using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
+using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel, Simplex_tree>;
 
 BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) {
   // ----------------------------------------------------------------------------
@@ -52,18 +45,6 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) {
   //
   // ----------------------------------------------------------------------------
   Point_cloud points;
-//   points.push_back({1., 0.});                  // 0
-//   points.push_back({0., 1.});                  // 1
-//   points.push_back({2., 1.});                  // 2
-//   points.push_back({3., 2.});                  // 3
-//   points.push_back({0., 3.});                  // 4
-//   points.push_back({3. + std::sqrt(3.), 3.});  // 5
-//   points.push_back({1., 4.});                  // 6
-//   points.push_back({3., 4.});                  // 7
-//   points.push_back({2., 4. + std::sqrt(3.)});  // 8
-//   points.push_back({0., 4.});                  // 9
-//   points.push_back({-0.5, 2.});                // 10
-
 
   std::vector<FT> point0({1., 0.});
   points.emplace_back(point0.begin(), point0.end());
@@ -156,36 +137,32 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) {
   for (std::size_t vertex = 0; vertex <= 2; vertex++) {
     points012.push_back(cech_complex_for_doc.get_point(vertex));
   }
-//   std::size_t dimension = points[0].end() - points[0].begin();
-//   Min_sphere ms012(dimension, points012.begin(), points012.end());
 
   Kernel kern;
   Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2}));
-  CGAL::NT_converter<FT, double> cast_to_double;
-  std::clog << "f012= " << f012 << " | points012_radius= " << std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points012.begin(), points012.end()))) << std::endl;
-
+  std::clog << "f012= " << f012 << std::endl;
 
+  CGAL::NT_converter<FT, double> cast_to_double;
   GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points012.begin(), points012.end()))));
 
   Point_cloud points1410;
   points1410.push_back(cech_complex_for_doc.get_point(1));
   points1410.push_back(cech_complex_for_doc.get_point(4));
   points1410.push_back(cech_complex_for_doc.get_point(10));
-//   Min_sphere ms1410(dimension, points1410.begin(), points1410.end());
 
   Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10}));
-  std::clog << "f1410= " << f1410 << " | points1410_radius= " << std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points1410.begin(), points1410.end()))) << std::endl;
+  std::clog << "f1410= " << f1410 << std::endl;
 
-  GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points1410.begin(), points1410.end()))));
+  // In this case, the computed sphere using CGAL kernel does not match the minimal enclosing ball; the filtration value check is therefore done against a hardcoded value
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, 1.);
 
   Point_cloud points469;
   points469.push_back(cech_complex_for_doc.get_point(4));
   points469.push_back(cech_complex_for_doc.get_point(6));
   points469.push_back(cech_complex_for_doc.get_point(9));
-//   Min_sphere ms469(dimension, points469.begin(), points469.end());
 
   Simplex_tree::Filtration_value f469 = st2.filtration(st2.find({4, 6, 9}));
-  std::clog << "f469= " << f469 << " | points469_radius= " << std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points469.begin(), points469.end()))) << std::endl;
+  std::clog << "f469= " << f469 << std::endl;
 
   GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points469.begin(), points469.end()))));
 
diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp
index ccd7d453..0c945cad 100644
--- a/src/Cech_complex/utilities/cech_persistence.cpp
+++ b/src/Cech_complex/utilities/cech_persistence.cpp
@@ -25,14 +25,12 @@
 // Types definition
 using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
 using Filtration_value = Simplex_tree::Filtration_value;
-// using Point = std::vector<double>;
-// using Point_cloud = std::vector<Point>;
 
 using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
 using Point = typename Kernel::Point_d;
 using Point_cloud = std::vector<Point>;
 using Points_off_reader = Gudhi::Points_off_reader<Point>;
-using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel>;
+using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel, Simplex_tree>;
 using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
 using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>;
 
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index f69ed6ec..85790baf 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -1278,10 +1278,8 @@ class Simplex_tree {
           intersection.emplace_back(next->first, Node(nullptr, filt));
         }
       }
-      std::clog << "Hind: after intersection insertion" << std::endl;
       if (intersection.size() != 0) {
         // Reverse the order to insert
-        std::clog << "Hind: declare new siblings" << std::endl;
         Siblings * new_sib = new Siblings(siblings,  // oncles
                                           simplex->first,  // parent
                                           boost::adaptors::reverse(intersection));  // boost::container::ordered_unique_range_t
@@ -1290,12 +1288,10 @@ class Simplex_tree {
         for (auto new_sib_member = new_sib->members().begin();
              new_sib_member != new_sib->members().end();
              new_sib_member++) {
-           std::clog << "Hind: check the blocker result" << std::endl;
            bool blocker_result = block_simplex(new_sib_member);
            // new_sib member has been blocked by the blocker function
            // add it to the list to be removed - do not perform it while looping on it
            if (blocker_result) {
-             std::clog << "Hind: add to list of blocked sib to be removed" << std::endl;
              blocked_new_sib_vertex_list.push_back(new_sib_member->first);
            }
         }
diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h
index ae5168aa..a8ee4a75 100644
--- a/src/common/include/gudhi/distance_functions.h
+++ b/src/common/include/gudhi/distance_functions.h
@@ -65,19 +65,17 @@ class Euclidean_distance {
  * The points are assumed to have the same dimension. */
 class Minimal_enclosing_ball_radius {
  public:
-   /** \brief TODO
+   /** \brief Enclosing ball radius from two points using CGAL.
    *
    * @param[in] point_1
    * @param[in] point_2
-   * @return
-   * \tparam Point
+   * @return Enclosing ball radius for the two points.
+   * \tparam Point must be a Kernel::Point_d from CGAL.
    *
    */
-  //typename FT = typename Kernel::FT,
   template< typename Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>,
             typename Point= typename Kernel::Point_d>
   double operator()(const Point& point_1, const Point& point_2) const {
-    std::clog << "Added template: distance betw points 1 and 2" << std::endl;
     Kernel kernel_;
     return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.;
   }
@@ -94,25 +92,21 @@ class Minimal_enclosing_ball_radius {
   template< typename Point >
   typename std::iterator_traits<typename boost::range_iterator<Point>::type>::value_type
   operator()(const Point& point_1, const Point& point_2) const {
-    std::clog << "Hind: Minimal_enclosing_ball_radius point1 et 2; Euclidean" << std::endl;
     std::clog << "#" << *point_1.begin() << "##" << *point_2.begin() << std::endl;
     return Euclidean_distance()(point_1, point_2) / 2.;
   }
 
-
-  /** \brief TODO
+  /** \brief Enclosing ball radius from a point cloud using CGAL.
    *
    * @param[in] point_cloud The points.
-   * @return
-   * \tparam Point_cloud
+   * @return Enclosing ball radius for the points.
+   * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL.
    *
    */
-  //typename FT = typename Kernel::FT,
   template< typename Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>,
             typename Point= typename Kernel::Point_d,
             typename Point_cloud = std::vector<Point>>
   double operator()(const Point_cloud& point_cloud) const {
-    std::clog << "Added template: distance in point cloud" << std::endl;
     Kernel kernel_;
     return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end())));
   }
@@ -133,8 +127,6 @@ class Minimal_enclosing_ball_radius {
             typename Coordinate = typename std::iterator_traits<Coordinate_iterator>::value_type>
   Coordinate
   operator()(const Point_cloud& point_cloud) const {
-    std::clog << "Hind: Minimal_enclosing_ball_radius point cloud; Miniball" << std::endl;
-
     using Min_sphere = Miniball::Miniball<Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
 
     Min_sphere ms(boost::size(*point_cloud.begin()), point_cloud.begin(), point_cloud.end());
diff --git a/src/common/include/gudhi/graph_simplicial_complex.h b/src/common/include/gudhi/graph_simplicial_complex.h
index 9190182c..da9dee7d 100644
--- a/src/common/include/gudhi/graph_simplicial_complex.h
+++ b/src/common/include/gudhi/graph_simplicial_complex.h
@@ -18,8 +18,6 @@
 #include <map>
 #include <tuple>  // for std::tie
 
-#include <iostream>
-
 namespace Gudhi {
 /** @file
  * @brief Graph simplicial complex methods
@@ -78,7 +76,6 @@ Proximity_graph<SimplicialComplexForProximityGraph> compute_proximity_graph(
   for (auto it_u = points.begin(); it_u != points.end(); ++it_u) {
     idx_v = idx_u + 1;
     for (auto it_v = it_u + 1; it_v != points.end(); ++it_v, ++idx_v) {
-      std::clog << "#idx_u" << idx_u << "#idx_v " << idx_v << std::endl;
       fil = distance(*it_u, *it_v);
       if (fil <= threshold) {
         edges.emplace_back(idx_u, idx_v);