From 6d1bfe69ea473088eb2ca3ede626a2accaa58205 Mon Sep 17 00:00:00 2001 From: fgodi Date: Tue, 29 Nov 2016 10:51:43 +0000 Subject: 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 --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 8 ++-- .../include/gudhi/Persistence_graph.h | 38 ++++++++++++------- src/Bottleneck_distance/test/bottleneck_chrono.cpp | 2 +- .../test/bottleneck_unit_test.cpp | 43 +++++++++++----------- 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 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::infinity(); std::vector 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((idmax - idmin) / alpha); + long step = static_cast((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 -#include +#include #include 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 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 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 u; std::vector v; + int alive_count; }; template 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::infinity() && it->second - it->first > e) + if(it->second == std::numeric_limits::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::infinity() && it->second - it->first > e) + if(it->second == std::numeric_limits::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 (u.size() + v.size()); } +inline bool Persistence_graph::alive_match() const{ + return alive_count==0; +} + inline std::vector Persistence_graph::sorted_distances() const { - // could be optimized - std::set sorted_distances; - sorted_distances.emplace(0.); - for (int u_point_index = 0; u_point_index < size(); ++u_point_index) + std::vector 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::infinity()); - return std::vector(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 millisecs_t; millisecs_t duration(std::chrono::duration_cast(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 unif1(0.,upper_bound); - std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); + std::uniform_real_distribution unif2(upper_bound/10000.,upper_bound/100.); std::default_random_engine re; std::vector< std::pair > 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.); } -- cgit v1.2.3