From b638f095644a1cda0fac00e1b11f19b6cec50473 Mon Sep 17 00:00:00 2001 From: fgodi Date: Fri, 16 Dec 2016 10:44:16 +0000 Subject: 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 --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 71 ++++++++++++++-------- .../include/gudhi/Persistence_graph.h | 8 +-- 2 files changed, 49 insertions(+), 30 deletions(-) (limited to 'src/Bottleneck_distance/include') 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 -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::infinity()) - return std::numeric_limits::infinity(); - std::vector 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((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 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((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 +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::infinity()) + return std::numeric_limits::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::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 Persistence_graph::sorted_distances() const { std::vector 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); } -- cgit v1.2.3