Copyrigth Sophia + tests

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

Former-commit-id: 353baa0c191df44a3bc420efd4424de2f5c6859e

7 files changed, 192 insertions, 105 deletions

diff --git a/src/Bottleneck/example/random_diagrams.cpp b/src/Bottleneck/example/random_diagrams.cpp
index 71f152a6..a9d0968d 100644
--- a/src/Bottleneck/example/random_diagrams.cpp
+++ b/src/Bottleneck/example/random_diagrams.cpp
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Francois Godi
  *
- *    Copyright (C) 2015  INRIA Saclay (France)
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
diff --git a/src/Bottleneck/include/gudhi/Graph_matching.h b/src/Bottleneck/include/gudhi/Graph_matching.h
index ea47e1d5..4cd3180e 100644
--- a/src/Bottleneck/include/gudhi/Graph_matching.h
+++ b/src/Bottleneck/include/gudhi/Graph_matching.h
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Francois Godi
  *
- *    Copyright (C) 2015  INRIA Saclay (France)
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
@@ -34,11 +34,11 @@ namespace Gudhi {
 namespace bottleneck {
 
 template<typename Persistence_diagram1, typename Persistence_diagram2>
-double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.);
+double bottleneck_distance(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e = 0.);
 
 class Graph_matching {
  public:
-  Graph_matching(const Persistence_diagrams_graph& g);
+  explicit Graph_matching(const Persistence_diagrams_graph& g);
   Graph_matching& operator=(const Graph_matching& m);
   bool perfect() const;
   bool multi_augment();
@@ -52,7 +52,7 @@ class Graph_matching {
 
   Layered_neighbors_finder* layering() const;
   bool augment(Layered_neighbors_finder* layered_nf, int u_start_index, int max_depth);
-  void update(std::deque<int>& path);
+  void update(std::deque<int>* path);
 };
 
 Graph_matching::Graph_matching(const Persistence_diagrams_graph& g)
@@ -128,36 +128,37 @@ Layered_neighbors_finder* Graph_matching::layering() const {
 }
 
 bool Graph_matching::augment(Layered_neighbors_finder *layered_nf, int u_start_index, int max_depth) {
-  std::deque<int> path;
-  path.emplace_back(u_start_index);
+  std::deque<int>* path = new std::deque<int>();
+  path->emplace_back(u_start_index);
   // start is a point from U
   do {
-    if (static_cast<int>(path.size()) > max_depth) {
-      path.pop_back();
-      path.pop_back();
+    if (static_cast<int>(path->size()) > max_depth) {
+      path->pop_back();
+      path->pop_back();
     }
-    if (path.empty())
+    if (path->empty())
       return false;
-    int w = path.back();
-    path.emplace_back(layered_nf->pull_near(w, path.size() / 2));
-    while (path.back() == null_point_index()) {
-      path.pop_back();
-      path.pop_back();
-      if (path.empty())
+    int w = path->back();
+    path->emplace_back(layered_nf->pull_near(w, path->size() / 2));
+    while (path->back() == null_point_index()) {
+      path->pop_back();
+      path->pop_back();
+      if (path->empty())
         return false;
-      path.pop_back();
-      path.emplace_back(layered_nf->pull_near(path.back(), path.size() / 2));
+      path->pop_back();
+      path->emplace_back(layered_nf->pull_near(path->back(), path->size() / 2));
     }
-    path.emplace_back(v_to_u.at(path.back()));
-  } while (path.back() != null_point_index());
-  path.pop_back();
+    path->emplace_back(v_to_u.at(path->back()));
+  } while (path->back() != null_point_index());
+  path->pop_back();
   update(path);
+  delete path;
   return true;
 }
 
-void Graph_matching::update(std::deque<int>& path) {
-  unmatched_in_u.remove(path.front());
-  for (auto it = path.cbegin(); it != path.cend(); ++it) {
+void Graph_matching::update(std::deque<int> *path) {
+  unmatched_in_u.remove(path->front());
+  for (auto it = path->cbegin(); it != path->cend(); ++it) {
     int tmp = *it;
     ++it;
     v_to_u[*it] = tmp;
@@ -165,14 +166,14 @@ void Graph_matching::update(std::deque<int>& path) {
 }
 
 template<typename Persistence_diagram1, typename Persistence_diagram2>
-double bottleneck_distance(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e) {
+double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) {
   Persistence_diagrams_graph g(diag1, diag2, e);
   std::vector<double>* sd = g.sorted_distances();
   int idmin = 0;
   int idmax = sd->size() - 1;
   double alpha = pow(sd->size(), 0.25);
   Graph_matching m(g);
-  Graph_matching biggest_unperfect = m;
+  Graph_matching biggest_unperfect(g);
   while (idmin != idmax) {
     int pas = static_cast<int>((idmax - idmin) / alpha);
     m.set_r(sd->at(idmin + pas));
diff --git a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
index de36e00b..80e5b5ed 100644
--- a/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
+++ b/src/Bottleneck/include/gudhi/Layered_neighbors_finder.h
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Francois Godi
  *
- *    Copyright (C) 2015  INRIA Saclay (France)
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
diff --git a/src/Bottleneck/include/gudhi/Neighbors_finder.h b/src/Bottleneck/include/gudhi/Neighbors_finder.h
index 98256571..6a05e1a0 100644
--- a/src/Bottleneck/include/gudhi/Neighbors_finder.h
+++ b/src/Bottleneck/include/gudhi/Neighbors_finder.h
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Francois Godi
  *
- *    Copyright (C) 2015  INRIA Saclay (France)
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
diff --git a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
index 7e278209..8e9eaaf6 100644
--- a/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
+++ b/src/Bottleneck/include/gudhi/Persistence_diagrams_graph.h
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Francois Godi
  *
- *    Copyright (C) 2015  INRIA Saclay (France)
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
@@ -44,100 +44,101 @@ int null_point_index();
 
 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).
-  template<typename Persistence_diagram1, typename Persistence_diagram2>
-  Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e = 0.);
-  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::vector<double>* sorted_distances();
+    // 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).
+    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::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;
+    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;
 };
 
 inline int null_point_index() {
-  return -1;
+    return -1;
 }
 
 template<typename Persistence_diagram1, typename Persistence_diagram2>
-Persistence_diagrams_graph::Persistence_diagrams_graph(Persistence_diagram1& diag1, Persistence_diagram2& diag2, double e)
+Persistence_diagrams_graph::Persistence_diagrams_graph(const Persistence_diagram1 &diag1,
+                                                       const Persistence_diagram2 &diag2, double e)
     : u(), v() {
-  for (auto it = diag1.cbegin(); it != diag1.cend(); ++it)
-    if (it->second - it->first > e)
-      u.emplace_back(*it);
-  for (auto it = diag2.cbegin(); it != diag2.cend(); ++it)
-    if (it->second - it->first > e)
-      v.emplace_back(*it);
-  if (u.size() < v.size())
-    swap(u, v);
+    for (auto it = diag1.cbegin(); it != diag1.cend(); ++it)
+        if (it->second - it->first > e)
+            u.emplace_back(*it);
+    for (auto it = diag2.cbegin(); it != diag2.cend(); ++it)
+        if (it->second - it->first > e)
+            v.emplace_back(*it);
+    if (u.size() < v.size())
+        swap(u, v);
 }
 
 Persistence_diagrams_graph::Persistence_diagrams_graph::Persistence_diagrams_graph()
     : u(), v() { }
 
 inline bool Persistence_diagrams_graph::on_the_u_diagonal(int u_point_index) const {
-  return u_point_index >= static_cast<int> (u.size());
+    return u_point_index >= static_cast<int> (u.size());
 }
 
 inline bool Persistence_diagrams_graph::on_the_v_diagonal(int v_point_index) const {
-  return v_point_index >= static_cast<int> (v.size());
+    return v_point_index >= static_cast<int> (v.size());
 }
 
 inline int Persistence_diagrams_graph::corresponding_point_in_u(int v_point_index) const {
-  return on_the_v_diagonal(v_point_index) ?
-      v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size());
+    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 {
-  return on_the_u_diagonal(u_point_index) ?
-      u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size());
+    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));
+    // 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 int Persistence_diagrams_graph::size() const {
-  return static_cast<int> (u.size() + v.size());
+    return static_cast<int> (u.size() + v.size());
 }
 
 inline 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)
-    for (int v_point_index = 0; v_point_index < size(); ++v_point_index)
-      sorted_distances.emplace(distance(u_point_index, v_point_index));
-  return new std::vector<double>(sorted_distances.cbegin(), sorted_distances.cend());
+    // could be optimized
+    std::set<double> sorted_distances;
+    for (int u_point_index = 0; u_point_index < size(); ++u_point_index)
+        for (int v_point_index = 0; v_point_index < size(); ++v_point_index)
+            sorted_distances.emplace(distance(u_point_index, v_point_index));
+    return new std::vector<double>(sorted_distances.cbegin(), sorted_distances.cend());
 }
 
 inline Diagram_point Persistence_diagrams_graph::get_u_point(int u_point_index) const {
-  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));
-  double x = (projector.first + projector.second) / 2;
-  return Diagram_point(x, x);
+    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));
+    double x = (projector.first + projector.second) / 2;
+    return Diagram_point(x, x);
 }
 
 inline Diagram_point Persistence_diagrams_graph::get_v_point(int v_point_index) const {
-  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));
-  double x = (projector.first + projector.second) / 2;
-  return Diagram_point(x, x);
+    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));
+    double x = (projector.first + projector.second) / 2;
+    return Diagram_point(x, x);
 }
 
 }  // namespace bottleneck
diff --git a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
index 4af672e4..3d04aa97 100644
--- a/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
+++ b/src/Bottleneck/include/gudhi/Planar_neighbors_finder.h
@@ -4,7 +4,7 @@
  *
  *    Author(s):       Francois Godi
  *
- *    Copyright (C) 2015  INRIA Saclay (France)
+ *    Copyright (C) 2015  INRIA Sophia-Antipolis (France)
  *
  *    This program is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
@@ -24,7 +24,6 @@
 #define SRC_BOTTLENECK_INCLUDE_GUDHI_PLANAR_NEIGHBORS_FINDER_H_
 
 #include <list>
-#include <iostream>
 #include <set>
 
 #include "Persistence_diagrams_graph.h"
diff --git a/src/Bottleneck/test/bottleneck_unit_test.cpp b/src/Bottleneck/test/bottleneck_unit_test.cpp
index 068b8690..ea7f51f2 100644
--- a/src/Bottleneck/test/bottleneck_unit_test.cpp
+++ b/src/Bottleneck/test/bottleneck_unit_test.cpp
@@ -1,26 +1,112 @@
 #define BOOST_TEST_MODULE bottleneck test
 
 #include <boost/test/included/unit_test.hpp>
-
+#include <random>
 #include "gudhi/Graph_matching.h"
-#include <iostream>
 
 using namespace Gudhi::bottleneck;
 
-BOOST_AUTO_TEST_CASE(random_diagrams) {
-  int n = 100;
-  // Random construction
-  std::vector< std::pair<double, double> > v1, v2;
-  for (int i = 0; i < n; i++) {
-    int a = rand() % n;
-    v1.emplace_back(a, a + rand() % (n - a));
-    int b = rand() % n;
-    v2.emplace_back(b, b + rand() % (n - b));
-  }
-  // v1 and v2 are persistence diagrams containing each 100 randoms points.
-  double b = bottleneck_distance(v1, v2, 0);
-  //
-  std::cout << b << std::endl;
-  const double EXPECTED_DISTANCE = 98.5;
-  BOOST_CHECK(b == EXPECTED_DISTANCE);
+Persistence_diagrams_graph* random_graph_generator(){
+    int n1 = 100;
+    int n2 = 120;
+    double upper_bound = 80;
+    // 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 new Persistence_diagrams_graph(v1, v2, 0.);
+}
+
+BOOST_AUTO_TEST_CASE(global){
+    int n = 100;
+    // Random construction
+    std::vector< std::pair<double, double> > v1, v2;
+    for (int i = 0; i < n; i++) {
+        int a = rand() % n;
+        v1.emplace_back(a, a + rand() % (n - a));
+        int b = rand() % n;
+        v2.emplace_back(b, b + rand() % (n - b));
+    }
+    //
+    BOOST_CHECK(bottleneck_distance(v1, v2, 1.) == 98);
 }
+
+
+BOOST_AUTO_TEST_CASE(persistence_diagrams_graph) {
+    Persistence_diagrams_graph* g = random_graph_generator();
+    std::vector<double>* d = g->sorted_distances();
+    //
+    BOOST_CHECK(!g->on_the_u_diagonal(99));
+    BOOST_CHECK(!g->on_the_u_diagonal(100));
+    BOOST_CHECK(!g->on_the_u_diagonal(119));
+    BOOST_CHECK(g->on_the_u_diagonal(120));
+    BOOST_CHECK(!g->on_the_v_diagonal(99));
+    BOOST_CHECK(g->on_the_v_diagonal(100));
+    BOOST_CHECK(g->on_the_v_diagonal(119));
+    BOOST_CHECK(g->on_the_v_diagonal(120));
+    //
+    BOOST_CHECK(g->corresponding_point_in_u(0)==120);
+    BOOST_CHECK(g->corresponding_point_in_u(100)==0);
+    BOOST_CHECK(g->corresponding_point_in_v(0)==100);
+    BOOST_CHECK(g->corresponding_point_in_v(120)==0);
+    //
+    BOOST_CHECK(g->size()==220);
+    //
+    BOOST_CHECK((int) d->size() <= 220*220 - 100*120 + 1); // could be more strict
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,0))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,99))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,100))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,119))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,120))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(0,219))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,0))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,99))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,100))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,119))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,120))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(100,219))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,0))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,99))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,100))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,119))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,120))==1);
+    BOOST_CHECK(std::count(d->begin(), d->end(), g->distance(219,219))==1);
+    //
+    delete g;
+    delete d;
+}
+
+BOOST_AUTO_TEST_CASE(planar_nf) {
+    Persistence_diagrams_graph* g = random_graph_generator();
+    Planar_neighbors_finder pnf = Planar_neighbors_finder(*g,1.);
+    for(int v_point_index=0; v_point_index<100; v_point_index+=2)
+        pnf.add(v_point_index);
+    pnf.remove(0);
+    pnf.remove(1);
+    //
+    BOOST_CHECK(!pnf.contains(0));
+    BOOST_CHECK(!pnf.contains(1));
+    BOOST_CHECK(pnf.contains(2));
+    BOOST_CHECK(!pnf.contains(3));
+    //
+    int v_point_index_1 = pnf.pull_near(120/2);
+    BOOST_CHECK((v_point_index_1 == -1) || (g->distance(120/2,v_point_index_1)<1.));
+    BOOST_CHECK(!pnf.contains(v_point_index_1));
+    int v_point_index_2 = pnf.pull_near(120/2);
+    BOOST_CHECK((v_point_index_2 == -1) || (g->distance(120/2,v_point_index_2)<1.));
+    BOOST_CHECK(!pnf.contains(v_point_index_2));
+    BOOST_CHECK((v_point_index_2 != -1) || (v_point_index_1 == -1));
+    pnf.add(v_point_index_1);
+    BOOST_CHECK(pnf.contains(v_point_index_1));
+}
+