separation of exactly and approximate

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

Former-commit-id: db94c84c29afbb94c0b0710a53e6ee417a5618df

3 files changed, 50 insertions, 31 deletions

diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h
index 24a31ac1..c34ea933 100644
--- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h
+++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h
@@ -30,43 +30,62 @@ namespace Gudhi {
 
 namespace persistence_diagram {
 
-/** \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
- */
-template<typename Persistence_diagram1, typename Persistence_diagram2>
-double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=0.) {
-    Persistence_graph g(diag1, diag2, e);
-    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.)
-        sd = g.sorted_distances();
-    long idmin = 0;
-    long idmax = e==0. ? sd.size() - 1 : g.diameter_bound()/e + 1;
-    double alpha = std::pow(g.size(), 1./5.);
+double bottleneck_distance_approx(Persistence_graph& g, double e) {
+    double b_lower_bound = 0.;
+    double b_upper_bound = g.diameter_bound();
+    const double alpha = std::pow(g.size(), 1./5.);
     Graph_matching m(g);
     Graph_matching biggest_unperfect(g);
-    while (idmin != idmax) {
-        long step = static_cast<long>((idmax - idmin - 1)/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)
+    while (b_upper_bound - b_lower_bound > 2*e) {
+        double step = b_lower_bound + (b_upper_bound - b_lower_bound)/alpha;
+        m.set_r(step);
+        while (m.multi_augment()); //compute a maximum matching (in the graph corresponding to the current r)
         if (m.perfect()) {
-            idmax = idmin + step;
             m = biggest_unperfect;
+            b_upper_bound = step;
         } else {
             biggest_unperfect = m;
-            idmin = idmin + step + 1;
+            b_lower_bound = step;
         }
     }
-    b = std::max(b, e == 0. ? sd.at(idmin) : e*(idmin));
-    return b;
+    return (b_lower_bound + b_upper_bound)/2.;
 }
 
+double bottleneck_distance_exact(Persistence_graph& g) {
+    std::vector<double> sd  = g.sorted_distances();
+    long lower_bound_i = 0;
+    long upper_bound_i = sd.size()-1;
+    const double alpha = std::pow(g.size(), 1./5.);
+    Graph_matching m(g);
+    Graph_matching biggest_unperfect(g);
+    while (lower_bound_i != upper_bound_i) {
+        long step = lower_bound_i + static_cast<long>((upper_bound_i - lower_bound_i - 1)/alpha);
+        m.set_r(sd.at(step));
+        while (m.multi_augment()); //compute a maximum matching (in the graph corresponding to the current r)
+        if (m.perfect()) {
+            m = biggest_unperfect;
+            upper_bound_i = step;
+        } else {
+            biggest_unperfect = m;
+            lower_bound_i = step + 1;
+        }
+    }
+    return sd.at(lower_bound_i);
+}
 
+/** \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
+ */
+template<typename Persistence_diagram1, typename Persistence_diagram2>
+double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=0.) {
+    Persistence_graph g(diag1, diag2, e);
+    if(g.bottleneck_alive() == std::numeric_limits<double>::infinity())
+        return std::numeric_limits<double>::infinity();
+    double b = (e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e));
+    return std::max(g.bottleneck_alive(),b);
+}
 
 }  // namespace persistence_diagram
 
diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
index 2ee55995..45a4d586 100644
--- a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
+++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h
@@ -96,7 +96,7 @@ Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1,
     std::sort(v_alive.begin(), v_alive.end());
     if(u_alive.size() != v_alive.size())
         b_alive = std::numeric_limits<double>::infinity();
-    else for(auto it_u=u_alive.cbegin(), it_v=v_alive.cbegin();it_u != u_alive.cend();++it_u, ++it_v)
+    else for(auto it_u=u_alive.cbegin(), it_v=v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v)
         b_alive = std::max(b_alive, std::fabs(*it_u - *it_v));
 }
 
@@ -136,7 +136,7 @@ inline double Persistence_graph::bottleneck_alive() const{
 
 inline std::vector<double> Persistence_graph::sorted_distances() const {
     std::vector<double> distances;
-    distances.push_back(0.);
+    distances.push_back(0.); //for empty diagrams
     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)
@@ -150,7 +150,7 @@ inline Internal_point Persistence_graph::get_u_point(int u_point_index) const {
     if (!on_the_u_diagonal(u_point_index))
         return u.at(u_point_index);
     Internal_point projector = v.at(corresponding_point_in_v(u_point_index));
-    double m = (projector.x() + projector.y()) / 2;
+    double m = (projector.x() + projector.y()) / 2.;
     return Internal_point(m,m,u_point_index);
 }
 
@@ -158,7 +158,7 @@ inline Internal_point Persistence_graph::get_v_point(int v_point_index) const {
     if (!on_the_v_diagonal(v_point_index))
         return v.at(v_point_index);
     Internal_point projector = u.at(corresponding_point_in_u(v_point_index));
-    double m = (projector.x() + projector.y()) / 2;
+    double m = (projector.x() + projector.y()) / 2.;
     return Internal_point(m,m,v_point_index);
 }
 
diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp
index dcf98460..b46a19e4 100644
--- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp
+++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp
@@ -51,7 +51,7 @@ int main(){
         if(i%3==0)
             v2.emplace_back(std::max(a,b),std::max(a,b)+y);
     }
-    double epsilon = 0.0000000000000001;
+    double epsilon =  std::numeric_limits<double>::min();
     std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
     double b = bottleneck_distance(v1,v2, epsilon);
     std::chrono::steady_clock::time_point end = std::chrono::steady_clock::now();