Merge of bottleneck-glisse branch to speed up bottleneck computation

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@2335 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 164b0e8de1d89e5b056dd01e63d5cb6dae6c67e1

2 files changed, 49 insertions, 37 deletions

diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h
index e1708c5b..f51e22e9 100644
--- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h
+++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h
@@ -26,7 +26,8 @@
 #include <gudhi/Neighbors_finder.h>
 
 #include <vector>
-#include <list>
+#include <unordered_set>
+#include <algorithm>
 
 namespace Gudhi {
 
@@ -40,8 +41,6 @@ class Graph_matching {
  public:
   /** \internal \brief Constructor constructing an empty matching. */
   explicit Graph_matching(Persistence_graph &g);
-  /** \internal \brief Copy operator. */
-  Graph_matching& operator=(const Graph_matching& m);
   /** \internal \brief Is the matching perfect ? */
   bool perfect() const;
   /** \internal \brief Augments the matching with a maximal set of edge-disjoint shortest augmenting paths. */
@@ -50,12 +49,12 @@ class Graph_matching {
   void set_r(double r);
 
  private:
-  Persistence_graph& g;
+  Persistence_graph* gp;
   double r;
   /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */
   std::vector<int> v_to_u;
   /** \internal \brief All the unmatched points in U. */
-  std::list<int> unmatched_in_u;
+  std::unordered_set<int> unmatched_in_u;
 
   /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */
   Layered_neighbors_finder layering() const;
@@ -66,17 +65,9 @@ class Graph_matching {
 };
 
 inline Graph_matching::Graph_matching(Persistence_graph& g)
-    : g(g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u() {
+    : gp(&g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u(g.size()) {
   for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index)
-    unmatched_in_u.emplace_back(u_point_index);
-}
-
-inline Graph_matching& Graph_matching::operator=(const Graph_matching& m) {
-  g = m.g;
-  r = m.r;
-  v_to_u = m.v_to_u;
-  unmatched_in_u = m.unmatched_in_u;
-  return *this;
+    unmatched_in_u.insert(u_point_index);
 }
 
 inline bool Graph_matching::perfect() const {
@@ -88,12 +79,12 @@ inline bool Graph_matching::multi_augment() {
     return false;
   Layered_neighbors_finder layered_nf(layering());
   int max_depth = layered_nf.vlayers_number()*2 - 1;
-  double rn = sqrt(g.size());
+  double rn = sqrt(gp->size());
   // verification of a necessary criterion in order to shortcut if possible
   if (max_depth < 0 || (unmatched_in_u.size() > rn && max_depth >= rn))
     return false;
   bool successful = false;
-  std::list<int> tries(unmatched_in_u);
+  std::vector<int> tries(unmatched_in_u.cbegin(), unmatched_in_u.cend());
   for (auto it = tries.cbegin(); it != tries.cend(); it++)
     // 'augment' has side-effects which have to be always executed, don't change order
     successful = augment(layered_nf, *it, max_depth) || successful;
@@ -133,12 +124,12 @@ inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u
 }
 
 inline Layered_neighbors_finder Graph_matching::layering() const {
-  std::list<int> u_vertices(unmatched_in_u);
-  std::list<int> v_vertices;
-  Neighbors_finder nf(g, r);
-  for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index)
+  std::vector<int> u_vertices(unmatched_in_u.cbegin(), unmatched_in_u.cend());
+  std::vector<int> v_vertices;
+  Neighbors_finder nf(*gp, r);
+  for (int v_point_index = 0; v_point_index < gp->size(); ++v_point_index)
     nf.add(v_point_index);
-  Layered_neighbors_finder layered_nf(g, r);
+  Layered_neighbors_finder layered_nf(*gp, r);
   for (int layer = 0; !u_vertices.empty(); layer++) {
     // one layer is one step in the BFS
     for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) {
@@ -166,7 +157,8 @@ inline Layered_neighbors_finder Graph_matching::layering() const {
 }
 
 inline void Graph_matching::update(std::vector<int>& path) {
-  unmatched_in_u.remove(path.front());
+  // Must return 1.
+  unmatched_in_u.erase(path.front());
   for (auto it = path.cbegin(); it != path.cend(); ++it) {
     // Be careful, the iterator is incremented twice each time
     int tmp = *it;
diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
index cd5486f8..bdc47578 100644
--- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
+++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h
@@ -25,9 +25,6 @@
 
 // Inclusion order is important for CGAL patch
 #include <CGAL/Kd_tree.h>
-#include <CGAL/Kd_tree_node.h>
-#include <CGAL/Orthogonal_k_neighbor_search.h>
-#include <CGAL/Weighted_Minkowski_distance.h>
 #include <CGAL/Search_traits.h>
 
 #include <gudhi/Persistence_graph.h>
@@ -40,6 +37,33 @@ namespace Gudhi {
 
 namespace persistence_diagram {
 
+/** \internal \brief Variant of CGAL::Fuzzy_iso_box to ensure that the box ic closed. It isn't clear how necessary that is.
+ */
+struct Square_query {
+  typedef CGAL::Dimension_tag<2> D;
+  typedef Internal_point Point_d;
+  typedef double FT;
+  bool contains(Point_d p) const {
+    return std::abs(p.x()-c.x())<=size && std::abs(p.y()-c.y())<=size;
+  }
+  bool inner_range_intersects(CGAL::Kd_tree_rectangle<FT,D> const&r) const {
+    return
+      r.max_coord(0) >= c.x() - size &&
+      r.min_coord(0) <= c.x() + size &&
+      r.max_coord(1) >= c.y() - size &&
+      r.min_coord(1) <= c.y() + size;
+  }
+  bool outer_range_contains(CGAL::Kd_tree_rectangle<FT,D> const&r) const {
+    return
+      r.min_coord(0) >= c.x() - size &&
+      r.max_coord(0) <= c.x() + size &&
+      r.min_coord(1) >= c.y() - size &&
+      r.max_coord(1) <= c.y() + size;
+  }
+  Point_d c;
+  FT size;
+};
+
 /** \internal \brief data structure used to find any point (including projections) in V near to a query point from U
  * (which can be a projection).
  *
@@ -51,9 +75,7 @@ namespace persistence_diagram {
 class Neighbors_finder {
   typedef CGAL::Dimension_tag<2> D;
   typedef CGAL::Search_traits<double, Internal_point, const double*, Construct_coord_iterator, D> Traits;
-  typedef CGAL::Weighted_Minkowski_distance<Traits> Distance;
-  typedef CGAL::Orthogonal_k_neighbor_search<Traits, Distance> K_neighbor_search;
-  typedef K_neighbor_search::Tree Kd_tree;
+  typedef CGAL::Kd_tree<Traits> Kd_tree;
 
  public:
   /** \internal \brief Constructor taking the near distance definition as parameter. */
@@ -123,15 +145,13 @@ inline int Neighbors_finder::pull_near(int u_point_index) {
   } else {
     // Is the query point near to a V point in the plane ?
     Internal_point u_point = g.get_u_point(u_point_index);
-    std::array<double, 2> w = {
-      {1., 1.}
-    };
-    K_neighbor_search search(kd_t, u_point, 1, 0., true, Distance(0, 2, w.begin(), w.end()));
-    auto it = search.begin();
-    if (it == search.end() || g.distance(u_point_index, it->first.point_index) > r)
+    auto neighbor = kd_t.search_any_point(Square_query{u_point, r});
+    if(!neighbor)
       return null_point_index();
-    tmp = it->first.point_index;
-    kd_t.remove(g.get_v_point(tmp));
+    tmp = neighbor->point_index;
+    auto point = g.get_v_point(tmp);
+    int idx = point.point_index;
+    kd_t.remove(point, [idx](Internal_point const&p){return p.point_index == idx;});
   }
   return tmp;
 }