infinite points separately computed

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

Former-commit-id: 5d1cd63a1b8ff26feffce86a4febe2f42d52340d

3 files changed, 33 insertions, 27 deletions

diff --git a/src/Bottleneck_distance/concept/Persistence_diagram.h b/src/Bottleneck_distance/concept/Persistence_diagram.h
index d55f5573..e4766628 100644
--- a/src/Bottleneck_distance/concept/Persistence_diagram.h
+++ b/src/Bottleneck_distance/concept/Persistence_diagram.h
@@ -27,24 +27,20 @@ namespace Gudhi {
 
 namespace bottleneck_distance {
 
-/** \brief Concept of Diagram_point. get<0>() must return the birth of the corresponding component and get<1>() its death.
- * If the component stay alive, get<1>() must return std::numeric_limits<double>::infinity().
+/** \brief Concept of Diagram_point. std::get<0>(point) must return the birth of the corresponding component and std::get<1>(point) its death.
+ * A valid implementation of this concept is std::pair<double,double>.
+ * Birth must be 0 or positive, death should be larger than birth, death can be std::numeric_limits<double>::infinity() for components which stay alive.
  *
  * \ingroup bottleneck_distance
  */
-struct Diagram_point{
-    double get<int>();
-};
+struct Diagram_point;
 
 /** \brief Concept of persistence diagram. It's a range of Diagram_point.
+ * std::begin(diagram) and std::end(diagram) must return corresponding iterators.
  *
  * \ingroup bottleneck_distance
  */
-struct Persistence_Diagram
-{
-    iterator<Diagram_point> begin();
-    iterator<Diagram_point> end();
-};
+struct Persistence_Diagram;
 
 }  // namespace bottleneck_distance
 
diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h
index 20225e80..4a69fb96 100644
--- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h
+++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h
@@ -24,12 +24,13 @@
 #define BOTTLENECK_H_
 
 #include <gudhi/Graph_matching.h>
+#include <cmath>
 
 namespace Gudhi {
 
 namespace bottleneck_distance {
 
-/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams.
+/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams (see Concepts).
  * If the last parameter e is not 0 (default value if not explicited), you get an additive e-approximation.
  *
  * \ingroup bottleneck_distance
@@ -37,7 +38,8 @@ 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())
+    double b = g.bottleneck_alive();
+    if(b == std::numeric_limits<double>::infinity())
         return std::numeric_limits<double>::infinity();
     std::vector<double> sd;
     if(e == 0.)
@@ -45,7 +47,7 @@ double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &di
     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);
+    double alpha = std::pow(idmax, 0.25);
     Graph_matching m(g);
     Graph_matching biggest_unperfect(g);
     while (idmin != idmax) {
@@ -61,7 +63,8 @@ double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &di
             idmin = idmin + step + 1;
         }
     }
-    return e == 0. ? sd.at(idmin) : e*(idmin);
+    b = std::max(b, e == 0. ? sd.at(idmin) : e*(idmin));
+    return b;
 }
 
 
diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
index d31a4826..8accb926 100644
--- a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
+++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
@@ -55,7 +55,7 @@ public:
     /** \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;
+    double bottleneck_alive() 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 for the diameter of the convex hull of all non infinite points */
@@ -68,26 +68,29 @@ public:
 private:
     std::vector<Internal_point> u;
     std::vector<Internal_point> v;
-    int alive_count;
+    std::vector<double> u_alive;
+    std::vector<double> v_alive;
 };
 
 template<typename Persistence_diagram1, typename Persistence_diagram2>
 Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1,
                                      const Persistence_diagram2 &diag2, double e)
-    : u(), v(), alive_count(0)
+    : u(), v(), u_alive(), v_alive()
 {
     for (auto it = std::begin(diag1); it != std::end(diag1); ++it){
         if(std::get<1>(*it) == std::numeric_limits<double>::infinity())
-            alive_count++;
-        if (std::get<1>(*it) - std::get<0>(*it) > e)
+            u_alive.push_back(std::get<0>(*it));
+        else if (std::get<1>(*it) - std::get<0>(*it) > e)
             u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size()));
     }
     for (auto it = std::begin(diag2); it != std::end(diag2); ++it){
         if(std::get<1>(*it) == std::numeric_limits<double>::infinity())
-            alive_count--;
-        if (std::get<1>(*it) - std::get<0>(*it) > e)
+            v_alive.push_back(std::get<0>(*it));
+        else if (std::get<1>(*it) - std::get<0>(*it) > e)
             v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size()));
     }
+    std::sort(u_alive.begin(), u_alive.end());
+    std::sort(v_alive.begin(), v_alive.end());
     if (u.size() < v.size())
         swap(u, v);
 }
@@ -115,8 +118,6 @@ inline double Persistence_graph::distance(int u_point_index, int 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);
-    if(p_u.y() == p_v.y())
-        return std::fabs(p_u.x() - p_v.x());
     return std::max(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y()));
 }
 
@@ -124,8 +125,13 @@ 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 double Persistence_graph::bottleneck_alive() const{
+    if(u_alive.size() != v_alive.size())
+        return std::numeric_limits<double>::infinity();
+    double max = 0.;
+    for(auto it_u=u_alive.cbegin(), it_v=v_alive.cbegin();it_u != u_alive.cend();++it_u, ++it_v)
+        max = std::max(max, std::fabs(*it_u - *it_v));
+    return max;
 }
 
 inline std::vector<double> Persistence_graph::sorted_distances() const {
@@ -159,12 +165,13 @@ inline Internal_point Persistence_graph::get_v_point(int v_point_index) const {
 inline double Persistence_graph::diameter_bound() const {
     double max = 0.;
     for(auto it = u.cbegin(); it != u.cend(); it++)
-        max = std::max(max,it->y() == std::numeric_limits<double>::infinity() ? it->x() : it->y());
+        max = std::max(max, it->y());
     for(auto it = v.cbegin(); it != v.cend(); it++)
-        max = std::max(max,it->y() == std::numeric_limits<double>::infinity() ? it->x() : it->y());
+        max = std::max(max, it->y());
     return max;
 }
 
+
 }  // namespace bottleneck_distance
 
 }  // namespace Gudhi