better gestion of infinite points

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

Former-commit-id: 28cca4bf35541a677fe08d5f0a8b89dfe9017962

4 files changed, 53 insertions, 38 deletions

diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h
index 59df655e..20225e80 100644
--- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h
+++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h
@@ -37,17 +37,19 @@ namespace bottleneck_distance {
 template<typename Persistence_diagram1, typename Persistence_diagram2>
 double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=0.) {
     Persistence_graph g(diag1, diag2, e);
+    if(!g.alive_match())
+        return std::numeric_limits<double>::infinity();
     std::vector<double> sd;
     if(e == 0.)
         sd = g.sorted_distances();
-    int idmin = 0;
-    int idmax = e==0. ? sd.size() - 1 : g.diameter_bound()/e + 1;
+    long idmin = 0;
+    long idmax = e==0. ? sd.size() - 1 : g.diameter_bound()/e + 1;
     // alpha can change the complexity
     double alpha = pow(idmax, 0.25);
     Graph_matching m(g);
     Graph_matching biggest_unperfect(g);
     while (idmin != idmax) {
-        int step = static_cast<int>((idmax - idmin) / alpha);
+        long step = static_cast<long>((idmax - idmin) / alpha);
         m.set_r(e == 0. ?  sd.at(idmin + step) : e*(idmin + step));
         while (m.multi_augment());
         //The above while compute a maximum matching (according to the r setted before)
diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
index 506dba52..ebd3adaf 100644
--- a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
+++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
@@ -24,7 +24,7 @@
 #define PERSISTENCE_GRAPH_H_
 
 #include <vector>
-#include <set>
+#include <algorithm>
 #include <gudhi/Internal_point.h>
 
 namespace Gudhi {
@@ -39,7 +39,7 @@ namespace bottleneck_distance {
  */
 class Persistence_graph {
 public:
-    /** \internal \brief Initializer taking 2 Point (concept) ranges as parameters. */
+    /** \internal \brief Constructor taking 2 Persistence_Diagrams (concept) as parameters. */
     template<typename Persistence_diagram1, typename Persistence_diagram2>
     Persistence_graph(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 ? */
@@ -54,9 +54,11 @@ public:
     double distance(int u_point_index, int v_point_index) const;
     /** \internal \brief Returns size = |U| = |V|. */
     int size() const;
+    /** \internal \brief Is there as many infinite points (alive components) in both diagrams ? */
+    bool alive_match() const;
     /** \internal \brief Returns the O(n^2) sorted distances between the points. */
     std::vector<double> sorted_distances() const;
-    /** \internal \brief Returns an upper bound of the diameter of the convex hull */
+    /** \internal \brief Returns an upper bound for the diameter of the convex hull of all non infinite points */
     double diameter_bound() const;
     /** \internal \brief Returns the corresponding internal point */
     Internal_point get_u_point(int u_point_index) const;
@@ -66,18 +68,23 @@ public:
 private:
     std::vector<Internal_point> u;
     std::vector<Internal_point> v;
+    int alive_count;
 };
 
 template<typename Persistence_diagram1, typename Persistence_diagram2>
 Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1,
                                      const Persistence_diagram2 &diag2, double e)
-    : u(), v()
+    : u(), v(), alive_count(0)
 {
     for (auto it = diag1.cbegin(); it != diag1.cend(); ++it)
-        if (it->second != std::numeric_limits<double>::infinity() && it->second - it->first > e)
+        if(it->second == std::numeric_limits<double>::infinity())
+            alive_count++;
+        else if (it->second - it->first > e)
             u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size()));
     for (auto it = diag2.cbegin(); it != diag2.cend(); ++it)
-        if (it->second != std::numeric_limits<double>::infinity() && it->second - it->first > e)
+        if(it->second == std::numeric_limits<double>::infinity())
+            alive_count--;
+        else if (it->second - it->first > e)
             v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size()));
     if (u.size() < v.size())
         swap(u, v);
@@ -113,15 +120,20 @@ inline int Persistence_graph::size() const {
     return static_cast<int> (u.size() + v.size());
 }
 
+inline bool Persistence_graph::alive_match() const{
+    return alive_count==0;
+}
+
 inline std::vector<double> Persistence_graph::sorted_distances() const {
-    // could be optimized
-    std::set<double> sorted_distances;
-    sorted_distances.emplace(0.);
-    for (int u_point_index = 0; u_point_index < size(); ++u_point_index)
+    std::vector<double> distances;
+    distances.push_back(0.);
+    for (int u_point_index = 0; u_point_index < size(); ++u_point_index){
+        distances.push_back(distance(u_point_index, corresponding_point_in_v(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));
-    sorted_distances.emplace(std::numeric_limits<double>::infinity());
-    return std::vector<double>(sorted_distances.begin(),sorted_distances.end());
+            distances.push_back(distance(u_point_index, v_point_index));
+    }
+    std::sort(distances.begin(), distances.end());
+    return distances;
 }
 
 inline Internal_point Persistence_graph::get_u_point(int u_point_index) const {
diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp
index 8e823897..8bd41104 100644
--- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp
+++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp
@@ -52,7 +52,7 @@ int main(){
             v2.emplace_back(std::max(a,b),std::max(a,b)+y);
     }
     std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
-    double b = compute(v1,v2, 0.0001);
+    double b = compute(v1,v2, 0.00001);
     std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();
     typedef std::chrono::duration<int,std::milli> millisecs_t;
     millisecs_t duration(std::chrono::duration_cast<millisecs_t>(end-start));
diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
index 84101274..2e67d436 100644
--- a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
+++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp
@@ -70,25 +70,25 @@ BOOST_AUTO_TEST_CASE(persistence_graph){
     //
     BOOST_CHECK(g.size()==(n1+n2));
     //
-    BOOST_CHECK((int) d.size() <= (n1+n2)*(n1+n2) - n1*n2 + 2);
-    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((int) d.size() == (n1+n2)*(n1+n2) + n1 + n2 + 1);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,0))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n1-1))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n1))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n2-1))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,n2))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(0,(n1+n2)-1))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,0))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n1-1))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n1))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n2-1))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,n2))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance(n1,(n1+n2)-1))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,0))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n1-1))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n1))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n2-1))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,n2))>0);
+    BOOST_CHECK(std::count(d.begin(), d.end(), g.distance((n1+n2)-1,(n1+n2)-1))>0);
 }
 
 BOOST_AUTO_TEST_CASE(neighbors_finder) {
@@ -146,7 +146,7 @@ BOOST_AUTO_TEST_CASE(graph_matching) {
 
 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.);
+    std::uniform_real_distribution<double> unif2(upper_bound/10000.,upper_bound/100.);
     std::default_random_engine re;
     std::vector< std::pair<double, double> > v1, v2;
     for (int i = 0; i < n1; i++) {
@@ -162,5 +162,6 @@ BOOST_AUTO_TEST_CASE(global){
             v2.emplace_back(std::max(a,b),std::max(a,b)+y);
     }
     BOOST_CHECK(compute(v1, v2) <= upper_bound/100.);
-    BOOST_CHECK(compute(v1, v2, upper_bound/1000.) <= upper_bound/100. + upper_bound/1000.);
+    BOOST_CHECK(compute(v1, v2, upper_bound/10000.) <= upper_bound/100. + upper_bound/10000.);
+    BOOST_CHECK(std::abs(compute(v1, v2) - compute(v1, v2, upper_bound/10000.)) <=  upper_bound/10000.);
 }