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

2 files changed, 36 insertions, 123 deletions

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