persistence_diagrams_graph static + others things

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

Former-commit-id: 2562cbd49534373dae6d6d23ec52ba5984080729

6 files changed, 332 insertions, 291 deletions

diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h
index 9f36e936..17412e6c 100644
--- a/src/Bottleneck/include/gudhi/Graph_matching.h
+++ b/src/Bottleneck/include/gudhi/Graph_matching.h
@@ -22,7 +22,7 @@
 
 #ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_
 #define SRC_BOTTLENECK_INCLUDE_GUDHI_GRAPH_MATCHING_H_
-//#define DEBUG
+
 #include <deque>
 #include <list>
 #include <vector>
@@ -33,31 +33,50 @@ namespace Gudhi {
 
 namespace bottleneck {
 
+/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams.
+ *
+ *
+ *
+ * \ingroup bottleneck_distance
+ */
 template<typename Persistence_diagram1, typename Persistence_diagram2>
 double bottleneck_distance(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e = 0.);
 
+/** \internal \brief Structure representing a graph matching. The graph is a Persistence_diagrams_graph.
+ *
+ * \ingroup bottleneck_distance
+ */
 class Graph_matching {
 public:
-    explicit Graph_matching(const Persistence_diagrams_graph& g);
+    /** \internal \brief Constructor constructing an empty matching. */
+    explicit Graph_matching();
+    /** \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. */
     bool multi_augment();
+    /** \internal \brief Sets the maximum length of the edges allowed to be added in the matching, 0 initially. */
     void set_r(double r);
 
 private:
-    const Persistence_diagrams_graph& g;
     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;
 
+    /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */
     std::unique_ptr<Layered_neighbors_finder> layering() const;
+    /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */
     bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth);
+    /** \internal \brief Update the matching with the simple augmenting path given as parameter. */
     void update(std::deque<int> & path);
 };
 
-Graph_matching::Graph_matching(const Persistence_diagrams_graph& g)
-    : g(g), r(0), v_to_u(g.size(), null_point_index()), unmatched_in_u() {
-    for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index)
+Graph_matching::Graph_matching()
+    : r(0.), v_to_u(G::size(), null_point_index()), unmatched_in_u() {
+    for (int u_point_index = 0; u_point_index < G::size(); ++u_point_index)
         unmatched_in_u.emplace_back(u_point_index);
 }
 
@@ -68,56 +87,36 @@ Graph_matching& Graph_matching::operator=(const Graph_matching& m) {
     return *this;
 }
 
-/* inline */ bool Graph_matching::perfect() const {
-#ifdef DEBUG
-    std::cout << "  perfect? unmatched_in_u.size = " << unmatched_in_u.size() << std::endl << std::flush;
-#endif
+inline bool Graph_matching::perfect() const {
     return unmatched_in_u.empty();
 }
 
-/* inline */ bool Graph_matching::multi_augment() {
+inline bool Graph_matching::multi_augment() {
     if (perfect())
         return false;
-#ifdef DEBUG
-    std::cout << "  multi augment" << std::endl << std::flush;
-#endif
     Layered_neighbors_finder layered_nf = *layering();
-    double rn = sqrt(g.size());
     int max_depth = layered_nf.vlayers_number()*2 - 1;
-#ifdef DEBUG
-    std::cout<< "    nb_max_layer = " << max_depth << std::endl << std::flush;
-#endif
-    // verification of a necessary criterion
+    double rn = sqrt(G::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);
-    for (auto it = tries.cbegin(); it != tries.cend(); it++){
-        const bool tmp = augment(layered_nf, *it, max_depth); //force augment evaluation even if successful is already true
-        successful = successful || tmp;
-    }
+    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;
     return successful;
 }
 
-/* inline */ void Graph_matching::set_r(double r) {
+inline void Graph_matching::set_r(double r) {
     this->r = r;
 }
 
-bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) {
-#ifdef DEBUG
-    std::cout << "    augment" << std::endl << std::flush;
-#endif
-#ifdef DEBUG
-    std::cout << "      u_start_index = " << u_start_index << std::endl << std::flush;
-#endif
+inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) {
+    //V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices.
     std::deque<int> path;
     path.emplace_back(u_start_index);
-    // u_start is a point from U
     do {
-#ifdef DEBUG
-        std::cout << "      do" << std::endl << std::flush;
-        std::cout << "        path.size = " << static_cast<int>(path.size()) << std::endl << std::flush;
-#endif
         if (static_cast<int>(path.size()) > max_depth) {
             path.pop_back();
             path.pop_back();
@@ -134,28 +133,22 @@ bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_
             path.emplace_back(layered_nf.pull_near(path.back(), path.size() / 2));
         }
         path.emplace_back(v_to_u.at(path.back()));
-#ifdef DEBUG
-        std::cout << "        v_to_u = " << path.back() << std::endl << std::flush;
-#endif
     } while (path.back() != null_point_index());
+    //if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex
     path.pop_back();
     update(path);
     return true;
 }
 
-std::unique_ptr<Layered_neighbors_finder> Graph_matching::layering() const {
-#ifdef DEBUG
-    std::cout << "    layering" << std::endl << std::flush;
-#endif
-    bool end = false;
-    int layer = 0;
+inline std::unique_ptr<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)
+    Neighbors_finder nf(r);
+    for (int v_point_index = 0; v_point_index < G::size(); ++v_point_index)
         nf.add(v_point_index);
-    std::unique_ptr<Layered_neighbors_finder> layered_nf(new Layered_neighbors_finder(g, r));
-    while (!u_vertices.empty()) {
+    std::unique_ptr<Layered_neighbors_finder> layered_nf(new Layered_neighbors_finder(r));
+    for(int layer = 0; !u_vertices.empty(); layer++) {
+        // one layer is one step in the BFS
         for (auto it = u_vertices.cbegin(); it != u_vertices.cend(); ++it) {
             std::unique_ptr< std::list<int> > u_succ = std::move(nf.pull_all_near(*it));
             for (auto it = u_succ->cbegin(); it != u_succ->cend(); ++it) {
@@ -163,27 +156,27 @@ std::unique_ptr<Layered_neighbors_finder> Graph_matching::layering() const {
                 v_vertices.emplace_back(*it);
             }
         }
+        // When the above for finishes, we have progress of one half-step (from U to V) in the BFS
         u_vertices.clear();
-        for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) {
+        bool end = false;
+        for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++)
             if (v_to_u.at(*it) == null_point_index())
+                // we stop when a nearest exposed V vertex (from U exposed vertices) has been found
                 end = true;
             else
                 u_vertices.emplace_back(v_to_u.at(*it));
-        }
+        // When the above for finishes, we have progress of one half-step (from V to U) in the BFS
         if (end)
             return layered_nf;
         v_vertices.clear();
-        layer++;
     }
     return layered_nf;
 }
 
-void Graph_matching::update(std::deque<int>& path) {
-#ifdef DEBUG
-    std::cout << "      update" << std::endl << std::flush;
-#endif
+inline void Graph_matching::update(std::deque<int>& path) {
     unmatched_in_u.remove(path.front());
     for (auto it = path.cbegin(); it != path.cend(); ++it) {
+        // Be careful, the iterator is incremented twice each time
         int tmp = *it;
         ++it;
         v_to_u[*it] = tmp;
@@ -192,27 +185,28 @@ void Graph_matching::update(std::deque<int>& path) {
 
 template<typename Persistence_diagram1, typename Persistence_diagram2>
 double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) {
-    Persistence_diagrams_graph g(diag1, diag2, e);
-    std::unique_ptr< std::vector<double> > sd = std::move(g.sorted_distances());
+    G::initialize(diag1, diag2, e);
+    std::unique_ptr< std::vector<double> > sd = std::move(G::sorted_distances());
     int idmin = 0;
     int idmax = sd->size() - 1;
+    // alpha can be modified, this will change the complexity
     double alpha = pow(sd->size(), 0.25);
-    Graph_matching m(g);
-    Graph_matching biggest_unperfect(g);
+    Graph_matching m;
+    Graph_matching biggest_unperfect;
     while (idmin != idmax) {
-        int pas = static_cast<int>((idmax - idmin) / alpha);
-        m.set_r(sd->at(idmin + pas));
+        int step = static_cast<int>((idmax - idmin) / alpha);
+        m.set_r(sd->at(idmin + step));
         while (m.multi_augment());
+        // The above while compute a maximum matching (according to the r setted before)
         if (m.perfect()) {
-            idmax = idmin + pas;
+            idmax = idmin + step;
             m = biggest_unperfect;
         } else {
             biggest_unperfect = m;
-            idmin = idmin + pas + 1;
+            idmin = idmin + step + 1;
         }
     }
-    double b = sd->at(idmin);
-    return b;
+    return sd->at(idmin);
 }
 
 }  // namespace bottleneck
diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
index 1cf83393..e6b7f30d 100644
--- a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
+++ b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
@@ -27,44 +27,50 @@
 
 #include "Neighbors_finder.h"
 
-// Layered_neighbors_finder is a data structure used to find if a query point from U has neighbors in V in a given
-// vlayer of the vlayered persistence diagrams graph. V's points have to be added manually using their index.
-// A neighbor returned is automatically removed.
-
 namespace Gudhi {
 
 namespace bottleneck {
 
+/** \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) in a layered graph layer given as parmeter.
+ *
+ * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically removed.
+ *
+ * \ingroup bottleneck_distance
+ */
 class Layered_neighbors_finder {
- public:
-  Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r);
-  void add(int v_point_index, int vlayer);
-  int pull_near(int u_point_index, int vlayer);
-  int vlayers_number() const;
+public:
+    /** \internal \brief Constructor taking the near distance definition as parameter. */
+    Layered_neighbors_finder(double r);
+    /** \internal \brief A point added will be possibly pulled. */
+    void add(int v_point_index, int vlayer);
+    /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */
+    int pull_near(int u_point_index, int vlayer);
+    /** \internal \brief Returns the number of layers. */
+    int vlayers_number() const;
 
- private:
-  const Persistence_diagrams_graph& g;
-  const double r;
-  std::vector<Neighbors_finder> neighbors_finder;
+private:
+    const double r;
+    std::vector<Neighbors_finder> neighbors_finder;
 };
 
-Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_diagrams_graph& g, double r) :
-    g(g), r(r), neighbors_finder() { }
+Layered_neighbors_finder::Layered_neighbors_finder(double r) :
+    r(r), neighbors_finder() { }
 
-/* inline */ void Layered_neighbors_finder::add(int v_point_index, int vlayer) {
-  for (int l = neighbors_finder.size(); l <= vlayer; l++)
-    neighbors_finder.emplace_back(Neighbors_finder(g, r));
-  neighbors_finder.at(vlayer).add(v_point_index);
+inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) {
+    for (int l = neighbors_finder.size(); l <= vlayer; l++)
+        neighbors_finder.emplace_back(Neighbors_finder(r));
+    neighbors_finder.at(vlayer).add(v_point_index);
 }
 
-/* inline */ int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) {
-  if (static_cast<int> (neighbors_finder.size()) <= vlayer)
-    return null_point_index();
-  return neighbors_finder.at(vlayer).pull_near(u_point_index);
+inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) {
+    if (static_cast<int> (neighbors_finder.size()) <= vlayer)
+        return null_point_index();
+    return neighbors_finder.at(vlayer).pull_near(u_point_index);
 }
 
-/* inline */ int Layered_neighbors_finder::vlayers_number() const {
-  return static_cast<int>(neighbors_finder.size());
+inline int Layered_neighbors_finder::vlayers_number() const {
+    return static_cast<int>(neighbors_finder.size());
 }
 
 }  // namespace bottleneck
diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h
index 0fad9889..ac3f8600 100644
--- a/src/Bottleneck/include/gudhi/Neighbors_finder.h
+++ b/src/Bottleneck/include/gudhi/Neighbors_finder.h
@@ -32,19 +32,25 @@ namespace Gudhi {
 
 namespace bottleneck {
 
-// Neighbors_finder is a data structure used to find if a query point from U has neighbors in V
-// in the persistence diagrams graph.
-// V's points have to be added manually using their index. A neighbor returned is automatically removed.
-
+/** \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).
+ *
+ * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically removed.
+ *
+ * \ingroup bottleneck_distance
+ */
 class Neighbors_finder {
 public:
-    Neighbors_finder(const Persistence_diagrams_graph& g, double r);
+    /** \internal \brief Constructor taking the near distance definition as parameter. */
+    Neighbors_finder(double r);
+    /** \internal \brief A point added will be possibly pulled. */
     void add(int v_point_index);
+    /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */
     int pull_near(int u_point_index);
+    /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */
     std::unique_ptr< std::list<int> > pull_all_near(int u_point_index);
 
 private:
-    const Persistence_diagrams_graph& g;
     const double r;
     Planar_neighbors_finder planar_neighbors_f;
     std::unordered_set<int> projections_f;
@@ -52,43 +58,46 @@ private:
     bool contains(int v_point_index);
 };
 
-Neighbors_finder::Neighbors_finder(const Persistence_diagrams_graph& g, double r) :
-    g(g), r(r), planar_neighbors_f(g, r), projections_f() { }
+Neighbors_finder::Neighbors_finder(double r) :
+    r(r), planar_neighbors_f(r), projections_f() { }
 
-/* inline */ void Neighbors_finder::add(int v_point_index) {
-    if (g.on_the_v_diagonal(v_point_index))
+inline void Neighbors_finder::add(int v_point_index) {
+    if (G::on_the_v_diagonal(v_point_index))
         projections_f.emplace(v_point_index);
     else
         planar_neighbors_f.add(v_point_index);
 }
 
-/* inline */ void Neighbors_finder::remove(int v_point_index) {
+inline void Neighbors_finder::remove(int v_point_index) {
     if(v_point_index == null_point_index())
         return;
-    if (g.on_the_v_diagonal(v_point_index))
+    if (G::on_the_v_diagonal(v_point_index))
         projections_f.erase(v_point_index);
     else
         planar_neighbors_f.remove(v_point_index);
 }
 
-/* inline */ bool Neighbors_finder::contains(int v_point_index) {
+inline bool Neighbors_finder::contains(int v_point_index) {
     return planar_neighbors_f.contains(v_point_index) || (projections_f.count(v_point_index)>0);
 }
 
-/* inline */ int Neighbors_finder::pull_near(int u_point_index) {
+inline int Neighbors_finder::pull_near(int u_point_index) {
     int tmp;
-    int c = g.corresponding_point_in_v(u_point_index);
-    if (g.on_the_u_diagonal(u_point_index) && projections_f.cbegin()!=projections_f.cend())
+    int c = G::corresponding_point_in_v(u_point_index);
+    if (G::on_the_u_diagonal(u_point_index) && !projections_f.empty())
+        //All projections are at distance 0
         tmp = *projections_f.cbegin();
-    else if (contains(c) && (g.distance(u_point_index, c) <= r))
+    else if (contains(c) && (G::distance(u_point_index, c) <= r))
+        //Is the query point near to its projection ?
         tmp = c;
     else
+        //Is the query point near to a V point in the plane ?
         tmp = planar_neighbors_f.pull_near(u_point_index);
     remove(tmp);
     return tmp;
 }
 
-std::unique_ptr< std::list<int> > Neighbors_finder::pull_all_near(int u_point_index) {
+inline std::unique_ptr< std::list<int> > Neighbors_finder::pull_all_near(int u_point_index) {
     std::unique_ptr< std::list<int> > all_pull = std::move(planar_neighbors_f.pull_all_near(u_point_index));
     int last_pull = pull_near(u_point_index);
     while (last_pull != null_point_index()) {
diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
index 644c1cd8..66e43145 100644
--- a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
+++ b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
@@ -26,52 +26,70 @@
 #include <vector>
 #include <set>
 #include <cmath>
-#include <utility>  // for pair<>
-#include <algorithm>  // for max
+#include <utility>
+#include <algorithm>
 
 namespace Gudhi {
 
 namespace bottleneck {
 
-// Diagram_point is the type of the persistence diagram's points
-typedef std::pair<double, double> Diagram_point;
-
-// Return the used index for encoding none of the points
+/** \internal \brief Returns the used index for encoding none of the points */
 int null_point_index();
 
-// Persistence_diagrams_graph is the interface beetwen any external representation of the two persistence diagrams and
-// the bottleneck distance computation. An interface is necessary to ensure basic functions complexity.
-
+/** \internal \brief Structure representing an euclidean bipartite graph containing
+ *  the points from the two persistence diagrams (including the projections).
+ *
+ * \ingroup bottleneck_distance
+ */
 class Persistence_diagrams_graph {
- public:
-    // Persistence_diagram1 and 2 are the types of any externals representations of persistence diagrams.
-    // They have to have an iterator over points, which have to have fields first (for birth) and second (for death).
+public:
+    /** \internal \brief Initializer taking 2 Point (concept) ranges as parameters. */
     template<typename Persistence_diagram1, typename Persistence_diagram2>
-    Persistence_diagrams_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e);
-    Persistence_diagrams_graph();
-    bool on_the_u_diagonal(int u_point_index) const;
-    bool on_the_v_diagonal(int v_point_index) const;
-    int corresponding_point_in_u(int v_point_index) const;
-    int corresponding_point_in_v(int u_point_index) const;
-    double distance(int u_point_index, int v_point_index) const;
-    int size() const;
-    std::unique_ptr< std::vector<double> > sorted_distances();
-
- private:
-    std::vector<Diagram_point> u;
-    std::vector<Diagram_point> v;
-    Diagram_point get_u_point(int u_point_index) const;
-    Diagram_point get_v_point(int v_point_index) const;
+    static void initialize(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e);
+    /** \internal \brief Is the given point from U the projection of a point in V ? */
+    static bool on_the_u_diagonal(int u_point_index);
+    /** \internal \brief Is the given point from V the projection of a point in U ? */
+    static bool on_the_v_diagonal(int v_point_index);
+    /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */
+    static int corresponding_point_in_u(int v_point_index);
+    /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */
+    static int corresponding_point_in_v(int u_point_index);
+    /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */
+    static double distance(int u_point_index, int v_point_index);
+    /** \internal \brief Returns size = |U| = |V|. */
+    static int size();
+    /** \internal \brief Returns the O(n^2) sorted distances between the points. */
+    static std::unique_ptr< std::vector<double> > sorted_distances();
+    /** \internal \brief Compare points regarding x coordinate. Use v_point_index for V points and u_point_index + G::size() for U points.  */
+    struct Compare_x{bool operator() (const int point_index_1, const int point_index_2) const;};
+    /** \internal \brief Compare points regarding y coordinate. Use v_point_index for V points and u_point_index + G::size() for U points. */
+    struct Compare_y{bool operator() (const int point_index_1, const int point_index_2) const;};
+
+private:
+    /** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */
+    typedef std::pair<double, double> Internal_point;
+    static std::vector<Internal_point> u;
+    static std::vector<Internal_point> v;
+    static Internal_point get_u_point(int u_point_index);
+    static Internal_point get_v_point(int v_point_index);
 };
 
-/* inline */ int null_point_index() {
+/** \internal \typedef \brief Alias used outside. */
+typedef Persistence_diagrams_graph G;
+
+// static initialization, seems to work but strange
+std::vector<G::Internal_point> Persistence_diagrams_graph::u = [] {return std::vector<G::Internal_point>();}();
+std::vector<G::Internal_point> Persistence_diagrams_graph::v = [] {return std::vector<G::Internal_point>();}();
+
+inline int null_point_index() {
     return -1;
 }
 
 template<typename Persistence_diagram1, typename Persistence_diagram2>
-Persistence_diagrams_graph::Persistence_diagrams_graph(const Persistence_diagram1 &diag1,
-                                                       const Persistence_diagram2 &diag2, double e)
-    : u(), v() {
+void Persistence_diagrams_graph::initialize(const Persistence_diagram1 &diag1,
+                                            const Persistence_diagram2 &diag2, double e){
+    u.clear();
+    v.clear();
     for (auto it = diag1.cbegin(); it != diag1.cend(); ++it)
         if (it->second - it->first > e)
             u.emplace_back(*it);
@@ -82,41 +100,37 @@ Persistence_diagrams_graph::Persistence_diagrams_graph(const Persistence_diagram
         swap(u, v);
 }
 
-Persistence_diagrams_graph::Persistence_diagrams_graph()
-    : u(), v() { }
-
-/* inline */ bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) const {
+inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) {
     return u_point_index >= static_cast<int> (u.size());
 }
 
-/* inline */ bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) const {
+inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) {
     return v_point_index >= static_cast<int> (v.size());
 }
 
-/* inline */ int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) const {
+inline int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) {
     return on_the_v_diagonal(v_point_index) ?
                 v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size());
 }
 
-/* inline */ int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) const {
+inline int Persistence_diagrams_graph::corresponding_point_in_v(int u_point_index) {
     return on_the_u_diagonal(u_point_index) ?
                 u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size());
 }
 
-/* inline */ double Persistence_diagrams_graph::distance(int u_point_index, int v_point_index) const {
-  // could be optimized for the case where one point is the projection of the other
-  if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index))
-    return 0;
-  Diagram_point p_u = get_u_point(u_point_index);
-  Diagram_point p_v = get_v_point(v_point_index);
-  return std::max(std::fabs(p_u.first - p_v.first), std::fabs(p_u.second - p_v.second));
+inline double Persistence_diagrams_graph::distance(int u_point_index, int v_point_index) {
+    if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index))
+        return 0;
+    Internal_point p_u = get_u_point(u_point_index);
+    Internal_point p_v = get_v_point(v_point_index);
+    return std::max(std::fabs(p_u.first - p_v.first), std::fabs(p_u.second - p_v.second));
 }
 
-/* inline */ int Persistence_diagrams_graph::size() const {
+inline int Persistence_diagrams_graph::size() {
     return static_cast<int> (u.size() + v.size());
 }
 
-/* inline */ std::unique_ptr< std::vector<double> > Persistence_diagrams_graph::sorted_distances() {
+inline std::unique_ptr< std::vector<double> > Persistence_diagrams_graph::sorted_distances() {
     // could be optimized
     std::set<double> sorted_distances;
     for (int u_point_index = 0; u_point_index < size(); ++u_point_index)
@@ -126,20 +140,32 @@ Persistence_diagrams_graph::Persistence_diagrams_graph()
     return sd_up;
 }
 
-/* inline */ Diagram_point Persistence_diagrams_graph::get_u_point(int u_point_index) const {
+inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::get_u_point(int u_point_index) {
     if (!on_the_u_diagonal(u_point_index))
         return u.at(u_point_index);
-    Diagram_point projector = v.at(corresponding_point_in_v(u_point_index));
+    Internal_point projector = v.at(corresponding_point_in_v(u_point_index));
     double x = (projector.first + projector.second) / 2;
-    return Diagram_point(x, x);
+    return Internal_point(x, x);
 }
 
-/* inline */ Diagram_point Persistence_diagrams_graph::get_v_point(int v_point_index) const {
+inline Persistence_diagrams_graph::Internal_point Persistence_diagrams_graph::get_v_point(int v_point_index) {
     if (!on_the_v_diagonal(v_point_index))
         return v.at(v_point_index);
-    Diagram_point projector = u.at(corresponding_point_in_u(v_point_index));
+    Internal_point projector = u.at(corresponding_point_in_u(v_point_index));
     double x = (projector.first + projector.second) / 2;
-    return Diagram_point(x, x);
+    return Internal_point(x, x);
+}
+
+inline bool Persistence_diagrams_graph::Compare_x::operator() (const int point_index_1, const int point_index_2) const{
+    G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size());
+    G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size());
+    return p1.first < p2.first;
+}
+
+inline bool Persistence_diagrams_graph::Compare_y::operator() (const int point_index_1, const int point_index_2) const{
+    G::Internal_point p1 = point_index_1 < G::size() ? G::get_v_point(point_index_1) : G::get_u_point(point_index_1 - G::size());
+    G::Internal_point p2 = point_index_2 < G::size() ? G::get_v_point(point_index_2) : G::get_u_point(point_index_2 - G::size());
+    return p1.second < p2.second;
 }
 
 }  // namespace bottleneck
diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
index e558239d..4d820c7f 100644
--- a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
+++ b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
@@ -27,88 +27,102 @@
 #include <set>
 
 #include "Persistence_diagrams_graph.h"
+#include "Grid_cell.h"
 
 namespace Gudhi {
 
 namespace bottleneck {
 
-// Planar_neighbors_finder is a data structure used to find if a query point from U has planar neighbors in V with the
-// planar distance.
-// V's points have to be added manually using their index. A neighbor returned is automatically removed but we can also
-// remove points manually using their index.
-
+/** \internal \brief Structure used to find any point in V near (according to the planar distance) to a query point from U.
+ *
+ * V points have to be added manually using their index and before the first remove/pull. A neighbor pulled is automatically removed. but we can also
+ * remove points manually using their index.
+ *
+ * \ingroup bottleneck_distance
+ */
 class Abstract_planar_neighbors_finder {
- public:
-  Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r);
-  virtual ~Abstract_planar_neighbors_finder() = 0;
-  virtual void add(int v_point_index) = 0;
-  virtual void remove(int v_point_index) = 0;
-  virtual bool contains(int v_point_index) const = 0;
-  virtual int pull_near(int u_point_index) = 0;
-  virtual std::unique_ptr< std::list<int> > pull_all_near(int u_point_index);
-
- protected:
-  const Persistence_diagrams_graph& g;
-  const double r;
+public:
+    /** \internal \brief Constructor TODO. */
+    Abstract_planar_neighbors_finder(double r);
+    virtual ~Abstract_planar_neighbors_finder() = 0;
+    /** \internal \brief A point added will be possibly pulled. */
+    virtual void add(int v_point_index) = 0;
+    /** \internal \brief A point manually removed will no longer be possibly pulled. */
+    virtual void remove(int v_point_index) = 0;
+    /** \internal \brief Can the point given as parameter be returned ? */
+    virtual bool contains(int v_point_index) const = 0;
+    /** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */
+    virtual int pull_near(int u_point_index) = 0;
+    /** \internal \brief Provide and remove all the V points near to the U point given as parameter. */
+    virtual std::unique_ptr< std::list<int> > pull_all_near(int u_point_index);
+
+protected:
+    const double r;
 };
 
-
-// Naive_pnf is a nave implementation of Abstract_planar_neighbors_finder
-
+/** \internal \brief Naive_pnf is an naïve Abstract_planar_neighbors_finder implementation
+ *
+ * \ingroup bottleneck_distance
+ */
 class Naive_pnf : public Abstract_planar_neighbors_finder {
- public:
-  Naive_pnf(const Persistence_diagrams_graph& g, double r);
-  void add(int v_point_index);
-  void remove(int v_point_index);
-  bool contains(int v_point_index) const;
-  int pull_near(int u_point_index);
-
- private:
-  std::set<int> candidates;
+public:
+    /** \internal \brief Constructor taking the near distance definition as parameter. */
+    Naive_pnf(double r);
+    /** \internal \brief A point added will be possibly pulled. */
+    void add(int v_point_index);
+    /** \internal \brief A point manually removed will no longer be possibly pulled. */
+    void remove(int v_point_index);
+    /** \internal \brief Can the point given as parameter be returned ? */
+    bool contains(int v_point_index) const;
+    /** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */
+    int pull_near(int u_point_index);
+
+private:
+    std::set<int> candidates;
 };
 
-
-// Planar_neighbors_finder is the used Abstract_planar_neighbors_finder's implementation
+/** \internal \typedef \brief Planar_neighbors_finder is the used Abstract_planar_neighbors_finder's implementation. */
 typedef Naive_pnf Planar_neighbors_finder;
 
-Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(const Persistence_diagrams_graph& g, double r) :
-    g(g), r(r) { }
 
-/* inline */ Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() { }
+Abstract_planar_neighbors_finder::Abstract_planar_neighbors_finder(double r) :
+    r(r) { }
 
-/* inline */ std::unique_ptr< std::list<int> > Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) {
-  std::unique_ptr< std::list<int> > all_pull(new std::list<int>);
-  int last_pull = pull_near(u_point_index);
-  while (last_pull != null_point_index()) {
-    all_pull->emplace_back(last_pull);
-    last_pull = pull_near(u_point_index);
-  }
-  return all_pull;
+inline Abstract_planar_neighbors_finder::~Abstract_planar_neighbors_finder() {}
+
+inline std::unique_ptr< std::list<int> > Abstract_planar_neighbors_finder::pull_all_near(int u_point_index) {
+    std::unique_ptr< std::list<int> > all_pull(new std::list<int>);
+    int last_pull = pull_near(u_point_index);
+    while (last_pull != null_point_index()) {
+        all_pull->emplace_back(last_pull);
+        last_pull = pull_near(u_point_index);
+    }
+    return all_pull;
 }
 
-Naive_pnf::Naive_pnf(const Persistence_diagrams_graph& g, double r) :
-    Abstract_planar_neighbors_finder(g, r), candidates() { }
+Naive_pnf::Naive_pnf(double r) :
+    Abstract_planar_neighbors_finder(r), candidates() { }
 
-/* inline */ void Naive_pnf::add(int v_point_index) {
-  candidates.emplace(v_point_index);
+inline void Naive_pnf::add(int v_point_index) {
+    candidates.emplace(v_point_index);
 }
 
-/* inline */ void Naive_pnf::remove(int v_point_index) {
-  candidates.erase(v_point_index);
+inline void Naive_pnf::remove(int v_point_index) {
+    candidates.erase(v_point_index);
 }
 
-/* inline */ bool Naive_pnf::contains(int v_point_index) const {
-  return (candidates.count(v_point_index) > 0);
+inline bool Naive_pnf::contains(int v_point_index) const {
+    return (candidates.count(v_point_index) > 0);
 }
 
-/* inline */ int Naive_pnf::pull_near(int u_point_index) {
-  for (auto it = candidates.begin(); it != candidates.end(); ++it)
-    if (g.distance(u_point_index, *it) <= r) {
-      int tmp = *it;
-      candidates.erase(it);
-      return tmp;
-    }
-  return null_point_index();
+inline int Naive_pnf::pull_near(int u_point_index) {
+    for (auto it = candidates.begin(); it != candidates.end(); ++it)
+        if (G::distance(u_point_index, *it) <= r) {
+            int tmp = *it;
+            candidates.erase(it);
+            return tmp;
+        }
+    return null_point_index();
 }
 
 }  // namespace bottleneck
diff --git a/src/Bottleneck/test/bottleneck_unit_test.cpp b/src/Bottleneck/test/bottleneck_unit_test.cpp
index 77f45867..926c8430 100644
--- a/src/Bottleneck/test/bottleneck_unit_test.cpp
+++ b/src/Bottleneck/test/bottleneck_unit_test.cpp
@@ -10,25 +10,6 @@ int n1 = 81; // a natural number >0
 int n2 = 180; // a natural number >0
 double upper_bound = 400.5; // any real >0
 
-std::unique_ptr<Persistence_diagrams_graph> random_graph_generator(){
-    // Random construction
-    std::uniform_real_distribution<double> unif(0.,upper_bound);
-    std::default_random_engine re;
-    std::vector< std::pair<double, double> > v1, v2;
-    for (int i = 0; i < n1; i++) {
-        double a = unif(re);
-        double b = unif(re);
-        v1.emplace_back(std::min(a,b), std::max(a,b));
-    }
-    for (int i = 0; i < n2; i++) {
-        double a = unif(re);
-        double b = unif(re);
-        v2.emplace_back(std::min(a,b), std::max(a,b));
-    }
-    return std::unique_ptr<Persistence_diagrams_graph>(new Persistence_diagrams_graph(v1, v2, 0.));
-}
-
-
 BOOST_AUTO_TEST_CASE(global){
     std::uniform_real_distribution<double> unif1(0.,upper_bound);
     std::uniform_real_distribution<double> unif2(upper_bound/1000.,upper_bound/100.);
@@ -49,51 +30,64 @@ BOOST_AUTO_TEST_CASE(global){
     BOOST_CHECK(bottleneck_distance(v1, v2) <= upper_bound/100.);
 }
 
-
-BOOST_AUTO_TEST_CASE(persistence_diagrams_graph) {
-    std::unique_ptr<Persistence_diagrams_graph> g = std::move(random_graph_generator());
-    std::unique_ptr< std::vector<double> > d = std::move(g->sorted_distances());
+BOOST_AUTO_TEST_CASE(persistence_diagrams_graph){
+    // Random construction
+    std::uniform_real_distribution<double> unif(0.,upper_bound);
+    std::default_random_engine re;
+    std::vector< std::pair<double, double> > v1, v2;
+    for (int i = 0; i < n1; i++) {
+        double a = unif(re);
+        double b = unif(re);
+        v1.emplace_back(std::min(a,b), std::max(a,b));
+    }
+    for (int i = 0; i < n2; i++) {
+        double a = unif(re);
+        double b = unif(re);
+        v2.emplace_back(std::min(a,b), std::max(a,b));
+    }
+    G::initialize(v1, v2, 0.);
+    std::unique_ptr< std::vector<double> > d = std::move(G::sorted_distances());
     //
-    BOOST_CHECK(!g->on_the_u_diagonal(n1-1));
-    BOOST_CHECK(!g->on_the_u_diagonal(n1));
-    BOOST_CHECK(!g->on_the_u_diagonal(n2-1));
-    BOOST_CHECK(g->on_the_u_diagonal(n2));
-    BOOST_CHECK(!g->on_the_v_diagonal(n1-1));
-    BOOST_CHECK(g->on_the_v_diagonal(n1));
-    BOOST_CHECK(g->on_the_v_diagonal(n2-1));
-    BOOST_CHECK(g->on_the_v_diagonal(n2));
+    BOOST_CHECK(!G::on_the_u_diagonal(n1-1));
+    BOOST_CHECK(!G::on_the_u_diagonal(n1));
+    BOOST_CHECK(!G::on_the_u_diagonal(n2-1));
+    BOOST_CHECK(G::on_the_u_diagonal(n2));
+    BOOST_CHECK(!G::on_the_v_diagonal(n1-1));
+    BOOST_CHECK(G::on_the_v_diagonal(n1));
+    BOOST_CHECK(G::on_the_v_diagonal(n2-1));
+    BOOST_CHECK(G::on_the_v_diagonal(n2));
     //
-    BOOST_CHECK(g->corresponding_point_in_u(0)==n2);
-    BOOST_CHECK(g->corresponding_point_in_u(n1)==0);
-    BOOST_CHECK(g->corresponding_point_in_v(0)==n1);
-    BOOST_CHECK(g->corresponding_point_in_v(n2)==0);
+    BOOST_CHECK(G::corresponding_point_in_u(0)==n2);
+    BOOST_CHECK(G::corresponding_point_in_u(n1)==0);
+    BOOST_CHECK(G::corresponding_point_in_v(0)==n1);
+    BOOST_CHECK(G::corresponding_point_in_v(n2)==0);
     //
-    BOOST_CHECK(g->size()==(n1+n2));
+    BOOST_CHECK(G::size()==(n1+n2));
     //
     BOOST_CHECK((int) d->size() <= (n1+n2)*(n1+n2) - n1*n2 + 1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,0))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n1-1))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n1))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n2-1))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,n2))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,(n1+n2)-1))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,0))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n1-1))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n1))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n2-1))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,n2))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(n1,(n1+n2)-1))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,0))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n1-1))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n1))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n2-1))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,n2))==1);
-    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance((n1+n2)-1,(n1+n2)-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,0))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n1-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n2-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,n2))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(0,(n1+n2)-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,0))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n1-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n2-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,n2))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance(n1,(n1+n2)-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,0))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n1-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n2-1))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,n2))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,(n1+n2)-1))==1);
 }
 
+
 BOOST_AUTO_TEST_CASE(planar_neighbors_finder) {
-    std::unique_ptr<Persistence_diagrams_graph> g = std::move(random_graph_generator());
-    Planar_neighbors_finder pnf = Planar_neighbors_finder(*g,1.);
+    Planar_neighbors_finder pnf = Planar_neighbors_finder(1.);
     for(int v_point_index=0; v_point_index<n1; v_point_index+=2)
         pnf.add(v_point_index);
     //
@@ -111,57 +105,51 @@ BOOST_AUTO_TEST_CASE(planar_neighbors_finder) {
     BOOST_CHECK(!pnf.contains(3));
     //
     int v_point_index_1 = pnf.pull_near(n2/2);
-    BOOST_CHECK((v_point_index_1 == -1) || (g->distance(n2/2,v_point_index_1)<=1.));
+    BOOST_CHECK((v_point_index_1 == -1) || (G::distance(n2/2,v_point_index_1)<=1.));
     BOOST_CHECK(!pnf.contains(v_point_index_1));
     std::list<int> l = *pnf.pull_all_near(n2/2);
     bool v = true;
     for(auto it = l.cbegin(); it != l.cend(); ++it)
-        v = v && (g->distance(n2/2,*it)>1.);
+        v = v && (G::distance(n2/2,*it)>1.);
     BOOST_CHECK(v);
     int v_point_index_2 = pnf.pull_near(n2/2);
     BOOST_CHECK(v_point_index_2 == -1);
-    pnf.add(v_point_index_1);
-    BOOST_CHECK(pnf.contains(v_point_index_1));
 }
 
 
 BOOST_AUTO_TEST_CASE(neighbors_finder) {
-    std::unique_ptr<Persistence_diagrams_graph> g = std::move(random_graph_generator());
-    Neighbors_finder nf = Neighbors_finder(*g,1.);
+    Neighbors_finder nf = Neighbors_finder(1.);
     for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2)
         nf.add(v_point_index);
     //
     int v_point_index_1 = nf.pull_near(n2/2);
-    BOOST_CHECK((v_point_index_1 == -1) || (g->distance(n2/2,v_point_index_1)<=1.));
+    BOOST_CHECK((v_point_index_1 == -1) || (G::distance(n2/2,v_point_index_1)<=1.));
     std::list<int> l = *nf.pull_all_near(n2/2);
     bool v = true;
     for(auto it = l.cbegin(); it != l.cend(); ++it)
-        v = v && (g->distance(n2/2,*it)>1.);
+        v = v && (G::distance(n2/2,*it)>1.);
     BOOST_CHECK(v);
     int v_point_index_2 = nf.pull_near(n2/2);
     BOOST_CHECK(v_point_index_2 == -1);
 }
 
 BOOST_AUTO_TEST_CASE(layered_neighbors_finder) {
-    std::unique_ptr<Persistence_diagrams_graph> g = std::move(random_graph_generator());
-    Layered_neighbors_finder lnf = Layered_neighbors_finder(*g,1.);
+    Layered_neighbors_finder lnf = Layered_neighbors_finder(1.);
     for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2)
         lnf.add(v_point_index, v_point_index % 7);
     //
     int v_point_index_1 = lnf.pull_near(n2/2,6);
-    BOOST_CHECK((v_point_index_1 == -1) || (g->distance(n2/2,v_point_index_1)<=1.));
+    BOOST_CHECK((v_point_index_1 == -1) || (G::distance(n2/2,v_point_index_1)<=1.));
     int v_point_index_2 = lnf.pull_near(n2/2,6);
     BOOST_CHECK(v_point_index_2 == -1);
     v_point_index_1 = lnf.pull_near(n2/2,0);
-    BOOST_CHECK((v_point_index_1 == -1) || (g->distance(n2/2,v_point_index_1)<=1.));
+    BOOST_CHECK((v_point_index_1 == -1) || (G::distance(n2/2,v_point_index_1)<=1.));
     v_point_index_2 = lnf.pull_near(n2/2,0);
     BOOST_CHECK(v_point_index_2 == -1);
 }
 
-
 BOOST_AUTO_TEST_CASE(graph_matching) {
-    std::unique_ptr<Persistence_diagrams_graph> g = std::move(random_graph_generator());
-    Graph_matching m1(*g);
+    Graph_matching m1;
     m1.set_r(0.);
     int e  = 0;
     while (m1.multi_augment())
@@ -177,3 +165,7 @@ BOOST_AUTO_TEST_CASE(graph_matching) {
     BOOST_CHECK(m2.perfect());
     BOOST_CHECK(!m1.perfect());
 }
+
+BOOST_AUTO_TEST_CASE(grid_cell) {
+
+}