diff options
6 files changed, 59 insertions, 143 deletions
diff --git a/src/Bipartite_graphs_matching/example/basic.cpp b/src/Bipartite_graphs_matching/example/basic.cpp index 26f9519e..772bcd4a 100644 --- a/src/Bipartite_graphs_matching/example/basic.cpp +++ b/src/Bipartite_graphs_matching/example/basic.cpp @@ -26,14 +26,14 @@ using namespace Gudhi::bipartite_graph_matching; int main() { - std::vector< Internal_point > v1, v2; + std::vector< std::pair<double,double> > v1, v2; - v1.push_back(Internal_point(2.7,3.7)); - v1.push_back(Internal_point(9.6,14)); - v1.push_back(Internal_point(34.2,34.974)); + v1.push_back(std::pair<double,double>(2.7,3.7)); + v1.push_back(std::pair<double,double>(9.6,14)); + v1.push_back(std::pair<double,double>(34.2,34.974)); - v2.push_back(Internal_point(2.8,4.45)); - v2.push_back(Internal_point(9.5,14.1)); + v2.push_back(std::pair<double,double>(2.8,4.45)); + v2.push_back(std::pair<double,double>(9.5,14.1)); double b = bottleneck_distance(v1, v2); diff --git a/src/Bipartite_graphs_matching/include/gudhi/Internal_point.h b/src/Bipartite_graphs_matching/include/gudhi/Internal_point.h index 633a3767..5d96ab93 100644 --- a/src/Bipartite_graphs_matching/include/gudhi/Internal_point.h +++ b/src/Bipartite_graphs_matching/include/gudhi/Internal_point.h @@ -23,9 +23,9 @@ #ifndef SRC_BOTTLENECK_INCLUDE_GUDHI_INTERNAL_POINT_H_ #define SRC_BOTTLENECK_INCLUDE_GUDHI_INTERNAL_POINT_H_ -//namespace Gudhi { +namespace Gudhi { -//namespace bipartite_graph_matching { +namespace bipartite_graph_matching { /** \internal \brief Returns the used index for encoding none of the points */ int null_point_index(); @@ -34,101 +34,25 @@ int null_point_index(); struct Internal_point { double vec[2]; int point_index; - Internal_point() { vec[0]= vec[1] = 0.; point_index = null_point_index(); } - Internal_point(double x, double y) { vec[0]=x; vec[1]=y; point_index = null_point_index(); } + Internal_point() {} + Internal_point(double x, double y, int p_i = null_point_index()) { vec[0]=x; vec[1]=y; point_index = p_i; } double x() const { return vec[ 0 ]; } double y() const { return vec[ 1 ]; } double& x() { return vec[ 0 ]; } double& y() { return vec[ 1 ]; } bool operator==(const Internal_point& p) const { - return (x() == p.x()) && (y() == p.y()); + return (x() == p.x()) && (y() == p.y()) && (point_index==p.point_index); } bool operator!=(const Internal_point& p) const { return ! (*this == p); } }; -namespace CGAL { -template <> -struct Kernel_traits<Internal_point> { - struct Kernel { - typedef double FT; - typedef double RT; - }; -}; -} - -struct Construct_coord_iterator { - typedef const double* result_type; - const double* operator()(const Internal_point& p) const - { return static_cast<const double*>(p.vec); } - const double* operator()(const Internal_point& p, int) const - { return static_cast<const double*>(p.vec+2); } -}; -/* -struct Distance { - typedef Internal_point Query_item; - typedef Internal_point Point_d; - typedef double FT; - typedef CGAL::Dimension_tag<2> D; - - FT transformed_distance(const Query_item& q, const Point_d& p) const { - FT d0= std::abs(q.x()-p.x()); - FT d1= std::abs(q.y()-p.y()); - return std::max(d0,d1); - } - FT min_distance_to_rectangle(const Query_item& q, - const CGAL::Kd_tree_rectangle<FT,D>& b) const { - FT d0(0.), d1(0.); - if (q.x() < b.min_coord(0)) - d0 += (b.min_coord(0)-q.x()); - if (q.x() > b.max_coord(0)) - d0 += (q.x()-b.max_coord(0)); - if (q.y() < b.min_coord(1)) - d1 += (b.min_coord(1)-q.y()); - if (q.y() > b.max_coord(1)) - d1 += (q.y()-b.max_coord(1)); - return std::max(d0,d1); - } - FT min_distance_to_rectangle(const Query_item& q, const CGAL::Kd_tree_rectangle<FT,D>& b,std::vector<FT>& dists){ - dists[0] = dists[1] = 0.; - if (q.x() < b.min_coord(0)) - dists[0] = (b.min_coord(0)- q.x()); - if (q.x() > b.max_coord(0)) - dists[0] = (q.x()-b.max_coord(0)); - if (q.y() < b.min_coord(1)) - dists[1] = (b.min_coord(1)-q.y()); - if (q.y() > b.max_coord(1)) - dists[1] = (q.y()-b.max_coord(1)); - return std::max(dists[0],dists[1]); - } - FT max_distance_to_rectangle(const Query_item& q, const CGAL::Kd_tree_rectangle<FT,D>& b) const { - FT d0 = (q.x() >= (b.min_coord(0)+b.max_coord(0))/2.) ? - (q.x()-b.min_coord(0)) : (b.max_coord(0)-q.x()); - FT d1 = (q.y() >= (b.min_coord(1)+b.max_coord(1))/2.) ? - (q.y()-b.min_coord(1)) : (b.max_coord(1)-q.y()); - return std::max(d0, d1); - } - FT max_distance_to_rectangle(const Query_item& q, const CGAL::Kd_tree_rectangle<FT,D>& b,std::vector<FT>& dists){ - dists[0] = (q.x() >= (b.min_coord(0)+b.max_coord(0))/2.0) ? - (q.x()-b.min_coord(0)) : (b.max_coord(0)-q.x()); - dists[1] = (q.y() >= (b.min_coord(1)+b.max_coord(1))/2.0) ? - (q.y()-b.min_coord(1)) : (b.max_coord(1)-q.y()); - return std::max(dists[0], dists[1]); - } - FT new_distance(FT& dist, FT old_off, FT new_off, - int) const { - return dist + new_off - old_off; //works ? - } - FT transformed_distance(FT d) const { return d; } - FT inverse_of_transformed_distance(FT d) { return d; } -}; // end of struct Distance -*/ inline int null_point_index() { return -1; } -//} // namespace bipartite_graph_matching +} // namespace bipartite_graph_matching -//} // namespace Gudhi +} // namespace Gudhi -#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_PERSISTENCE_DIAGRAMS_GRAPH_H_ +#endif // SRC_BOTTLENECK_INCLUDE_GUDHI_INTERNAL_POINT_H_ diff --git a/src/Bipartite_graphs_matching/include/gudhi/Neighbors_finder.h b/src/Bipartite_graphs_matching/include/gudhi/Neighbors_finder.h index 2ef7f768..bb6f5ea0 100644 --- a/src/Bipartite_graphs_matching/include/gudhi/Neighbors_finder.h +++ b/src/Bipartite_graphs_matching/include/gudhi/Neighbors_finder.h @@ -24,7 +24,6 @@ #define SRC_BOTTLENECK_INCLUDE_GUDHI_NEIGHBORS_FINDER_H_ #include <unordered_set> - #include "Planar_neighbors_finder.h" namespace Gudhi { diff --git a/src/Bipartite_graphs_matching/include/gudhi/Persistence_diagrams_graph.h b/src/Bipartite_graphs_matching/include/gudhi/Persistence_diagrams_graph.h index 2ba960a7..a460d814 100644 --- a/src/Bipartite_graphs_matching/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bipartite_graphs_matching/include/gudhi/Persistence_diagrams_graph.h @@ -85,11 +85,11 @@ inline void G::initialize(const Persistence_diagram1 &diag1, u.clear(); v.clear(); for (auto it = diag1.cbegin(); it != diag1.cend(); ++it) - if (it->y() - it->x() > e) - u.emplace_back(*it); + if (it->second - it->first > e) + u.push_back(Internal_point(it->first, it->second, u.size())); for (auto it = diag2.cbegin(); it != diag2.cend(); ++it) - if (it->y() - it->x() > e) - v.emplace_back(*it); + if (it->second - it->first > e) + v.push_back(Internal_point(it->first, it->second, v.size())); if (u.size() < v.size()) swap(u, v); } @@ -139,7 +139,7 @@ inline Internal_point G::get_u_point(int 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; - return Internal_point(m, m); + return Internal_point(m,m,u_point_index); } inline Internal_point G::get_v_point(int v_point_index) { @@ -147,7 +147,7 @@ inline Internal_point G::get_v_point(int 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; - return Internal_point(m, m); + return Internal_point(m,m,v_point_index); } } // namespace bipartite_graph_matching diff --git a/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h b/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h index 47c1a8fc..abb14e29 100644 --- a/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h +++ b/src/Bipartite_graphs_matching/include/gudhi/Planar_neighbors_finder.h @@ -28,6 +28,7 @@ #include <CGAL/Search_traits.h> #include <CGAL/Orthogonal_incremental_neighbor_search.h> #include <CGAL/Weighted_Minkowski_distance.h> +#include "../CGAL/Miscellaneous.h" #include "Persistence_diagrams_graph.h" @@ -66,9 +67,9 @@ private: class Cgal_pnf { typedef CGAL::Dimension_tag<2> D; - typedef CGAL::Search_traits<double, Internal_point, const double*, Construct_coord_iterator, D> Cgal_traits; - typedef CGAL::Weighted_Minkowski_distance<Cgal_traits> Distance; - typedef CGAL::Orthogonal_incremental_neighbor_search<Cgal_traits, Distance> K_neighbor_search; + typedef CGAL::Search_traits<double, Internal_point, const double*, CGAL::Construct_coord_iterator, D> Traits; + typedef CGAL::Weighted_Minkowski_distance<Traits> Distance; + typedef CGAL::Orthogonal_incremental_neighbor_search<Traits, Distance> K_neighbor_search; typedef K_neighbor_search::Tree Kd_tree; @@ -164,15 +165,14 @@ inline Cgal_pnf::Cgal_pnf(double r_) /** \internal \brief A point added will be possibly pulled. */ inline void Cgal_pnf::add(int v_point_index){ - Internal_point v_point = G::get_v_point(v_point_index); - v_point.point_index = v_point_index; - kd_t.insert(v_point); contents.insert(v_point_index); + kd_t.insert(G::get_v_point(v_point_index)); } /** \internal \brief A point manually removed will no longer be possibly pulled. */ inline void Cgal_pnf::remove(int v_point_index){ contents.erase(v_point_index); + kd_t.remove(G::get_v_point(v_point_index)); } /** \internal \brief Can the point given as parameter be returned ? */ @@ -183,23 +183,16 @@ inline bool Cgal_pnf::contains(int v_point_index) const{ /** \internal \brief Provide and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ inline int Cgal_pnf::pull_near(int u_point_index){ Internal_point u_point = G::get_u_point(u_point_index); - K_neighbor_search search(kd_t, u_point, 0., true, Distance(0.)); - for (auto it = search.begin(); it != search.end(); it++) - if(contents.count(it->first.point_index)==0) - kd_t.remove(it->first); - else if(G::distance(u_point_index, it->first.point_index) > r){ - for(auto itc=contents.cbegin(); itc != contents.cend(); itc++) - if(G::distance(u_point_index, *itc) <= r) - std::cout << G::distance(u_point_index, *itc) << " ! > " << r << std::endl; - return null_point_index(); - } - else - { - kd_t.remove(it->first); - contents.erase(it->first.point_index); - return it->first.point_index; - } - return null_point_index(); + std::vector<double> w = {1., 1.}; + K_neighbor_search search(kd_t, u_point, 0., true, Distance(0, 2, w)); + auto it = search.begin(); + if(it==search.end() || G::distance(u_point_index, it->first.point_index) > r){ + for(auto itc=contents.cbegin(); itc != contents.cend(); itc++) + if(G::distance(u_point_index, *itc) <= r) + std::cout << G::distance(u_point_index, *itc) << " ! > " << r << std::endl; + return null_point_index(); + } + return it->first.point_index; } inline std::shared_ptr< std::list<int> > Cgal_pnf::pull_all_near(int u_point_index) { diff --git a/src/Bipartite_graphs_matching/test/bottleneck_unit_test.cpp b/src/Bipartite_graphs_matching/test/bottleneck_unit_test.cpp index 2e79ccdc..d53a5a3c 100644 --- a/src/Bipartite_graphs_matching/test/bottleneck_unit_test.cpp +++ b/src/Bipartite_graphs_matching/test/bottleneck_unit_test.cpp @@ -13,31 +13,11 @@ int n1 = 81; // a natural number >0 int n2 = 180; // a natural number >0 double upper_bound = 406.43; // any real >0 -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::default_random_engine re; - std::vector< Internal_point > v1, v2; - for (int i = 0; i < n1; i++) { - double a = unif1(re); - double b = unif1(re); - double x = unif2(re); - double y = unif2(re); - v1.emplace_back(std::min(a,b), std::max(a,b)); - v2.emplace_back(std::min(a,b)+std::min(x,y), std::max(a,b)+std::max(x,y)); - if(i%5==0) - v1.emplace_back(std::min(a,b),std::min(a,b)+x); - if(i%3==0) - v2.emplace_back(std::max(a,b),std::max(a,b)+y); - } - BOOST_CHECK(bottleneck_distance(v1, v2) <= upper_bound/100.); -} - BOOST_AUTO_TEST_CASE(persistence_diagrams_graph){ // Random construction std::uniform_real_distribution<double> unif(0.,upper_bound); std::default_random_engine re; - std::vector< Internal_point > v1, v2; + std::vector< std::pair<double, double> > v1, v2; for (int i = 0; i < n1; i++) { double a = unif(re); double b = unif(re); @@ -88,7 +68,6 @@ BOOST_AUTO_TEST_CASE(persistence_diagrams_graph){ BOOST_CHECK(std::count(d->begin(), d->end(), G::distance((n1+n2)-1,(n1+n2)-1))==1); } - BOOST_AUTO_TEST_CASE(planar_neighbors_finder) { Planar_neighbors_finder pnf(1.); for(int v_point_index=0; v_point_index<n1; v_point_index+=2) @@ -167,6 +146,27 @@ BOOST_AUTO_TEST_CASE(graph_matching) { BOOST_CHECK(m2.perfect()); BOOST_CHECK(!m1.perfect()); } + +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::default_random_engine re; + std::vector< std::pair<double, double> > v1, v2; + for (int i = 0; i < n1; i++) { + double a = unif1(re); + double b = unif1(re); + double x = unif2(re); + double y = unif2(re); + v1.emplace_back(std::min(a,b), std::max(a,b)); + v2.emplace_back(std::min(a,b)+std::min(x,y), std::max(a,b)+std::max(x,y)); + if(i%5==0) + v1.emplace_back(std::min(a,b),std::min(a,b)+x); + if(i%3==0) + v2.emplace_back(std::max(a,b),std::max(a,b)+y); + } + BOOST_CHECK(bottleneck_distance(v1, v2) <= upper_bound/100.); +} + /* BOOST_AUTO_TEST_CASE(chrono) { std::ofstream objetfichier; @@ -176,7 +176,7 @@ BOOST_AUTO_TEST_CASE(chrono) { std::uniform_real_distribution<double> unif1(0.,upper_bound); std::uniform_real_distribution<double> unif2(upper_bound/1000.,upper_bound/100.); std::default_random_engine re; - std::vector< Internal_point > v1, v2; + std::vector< std::pair<double, double> > v1, v2; for (int i = 0; i < n; i++) { double a = unif1(re); double b = unif1(re); |