From 6dd7403da12dff125e1dedacbc2ed0ec67e643ba Mon Sep 17 00:00:00 2001 From: fgodi Date: Wed, 23 Nov 2016 14:55:25 +0000 Subject: Bottleneck.h git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1774 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 3b2f26c67ceb8a61f3efe3952f0618fff98465e4 --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 125 +++++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 src/Bottleneck_distance/include/gudhi/Bottleneck.h (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h new file mode 100644 index 00000000..6b6b1552 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -0,0 +1,125 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Francois Godi + * + * Copyright (C) 2015 INRIA Sophia-Antipolis (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef BOTTLENECK_H_ +#define BOTTLENECK_H_ + +#include + +namespace Gudhi { + +namespace bottleneck_distance { + +/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams. You get an additive e-approximation. + * + * + * \ingroup bottleneck_distance + */ +template +double compute(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e = 0.); + +template +double compute_exactly(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2); + +template +double compute_exactly(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2) { + G::initialize(diag1, diag2, 0.); + std::shared_ptr< std::vector > sd(G::sorted_distances()); + int idmin = 0; + int idmax = sd->size() - 1; + // alpha can be modified, this will change the complexity + double alpha = pow(sd->size(), 0.25); + Graph_matching m; + Graph_matching biggest_unperfect; + while (idmin != idmax) { + int step = static_cast((idmax - idmin) / alpha); + m.set_r(sd->at(idmin + step)); + while (m.multi_augment()); + //The above while compute a maximum matching (according to the r setted before) + if (m.perfect()) { + idmax = idmin + step; + m = biggest_unperfect; + } else { + biggest_unperfect = m; + idmin = idmin + step + 1; + } + } + return sd->at(idmin); +} + +template +double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { + if(e< std::numeric_limits::min()) + return compute_exactly(diag1, diag2); + G::initialize(diag1, diag2, e); + int in = G::diameter()/e; + int idmin = 0; + int idmax = in; + // alpha can be modified, this will change the complexity + double alpha = pow(in, 0.25); + Graph_matching m; + Graph_matching biggest_unperfect; + while (idmin != idmax) { + int step = static_cast((idmax - idmin) / alpha); + m.set_r(e*(idmin + step)); + while (m.multi_augment()); + //The above while compute a maximum matching (according to the r setted before) + if (m.perfect()) { + idmax = idmin + step; + m = biggest_unperfect; + } else { + biggest_unperfect = m; + idmin = idmin + step + 1; + } + } + return e*(idmin); +} + + + +} // namespace bottleneck_distance + +} // namespace Gudhi + +#endif // BOTTLENECK_H_ + +/* Dichotomic version +template +double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { + if(e< std::numeric_limits::min()) + return compute_exactly(diag1, diag2); + G::initialize(diag1, diag2, e); + double d = 0.; + double f = G::diameter(); + while (f-d > e){ + Graph_matching m; + m.set_r((d+f)/2.); + while (m.multi_augment()); + //The above while compute a maximum matching (according to the r setted before) + if (m.perfect()) + f = (d+f)/2.; + else + d= (d+f)/2.; + } + return d; +} */ + -- cgit v1.2.3 From 33ce062bf29b9bb7db9e83cde39cab6411d613dc Mon Sep 17 00:00:00 2001 From: fgodi Date: Wed, 23 Nov 2016 16:02:39 +0000 Subject: functions no longer return smart pointers git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1779 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: bc5e1b480774441c9af5e0ed902c455c6453c480 --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 10 +++++----- .../include/gudhi/CGAL/Kd_tree.h | 1 - .../include/gudhi/CGAL/Kd_tree_node.h | 3 ++- .../include/gudhi/Graph_matching.h | 22 ++++++++++------------ .../include/gudhi/Neighbors_finder.h | 8 ++++---- .../include/gudhi/Persistence_diagrams_graph.h | 12 +++--------- .../include/gudhi/Planar_neighbors_finder.h | 19 ++++++++----------- 7 files changed, 32 insertions(+), 43 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index 6b6b1552..1ae7788c 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -43,16 +43,16 @@ double compute_exactly(const Persistence_diagram1& diag1, const Persistence_diag template double compute_exactly(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2) { G::initialize(diag1, diag2, 0.); - std::shared_ptr< std::vector > sd(G::sorted_distances()); + std::vector sd(G::sorted_distances()); int idmin = 0; - int idmax = sd->size() - 1; + int idmax = sd.size() - 1; // alpha can be modified, this will change the complexity - double alpha = pow(sd->size(), 0.25); + double alpha = pow(sd.size(), 0.25); Graph_matching m; Graph_matching biggest_unperfect; while (idmin != idmax) { int step = static_cast((idmax - idmin) / alpha); - m.set_r(sd->at(idmin + step)); + m.set_r(sd.at(idmin + step)); while (m.multi_augment()); //The above while compute a maximum matching (according to the r setted before) if (m.perfect()) { @@ -63,7 +63,7 @@ double compute_exactly(const Persistence_diagram1 &diag1, const Persistence_diag idmin = idmin + step + 1; } } - return sd->at(idmin); + return sd.at(idmin); } template diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h index 7b1676cf..9db9b8b6 100644 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h +++ b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h @@ -28,7 +28,6 @@ #include #include -#include #include #include diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h index acefe926..d6d01439 100644 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h +++ b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h @@ -21,7 +21,8 @@ #ifndef CGAL_KD_TREE_NODE_H #define CGAL_KD_TREE_NODE_H -#include +#include "Splitters.h" + #include #include diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h index a8d68a9d..c50c3a9e 100644 --- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -23,8 +23,6 @@ #ifndef GRAPH_MATCHING_H_ #define GRAPH_MATCHING_H_ -#include - #include namespace Gudhi { @@ -56,11 +54,11 @@ private: std::list unmatched_in_u; /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ - std::shared_ptr layering() const; + Layered_neighbors_finder layering() const; /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */ bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth); /** \internal \brief Update the matching with the simple augmenting path given as parameter. */ - void update(std::deque & path); + void update(std::vector & path); }; inline Graph_matching::Graph_matching() @@ -83,7 +81,7 @@ inline bool Graph_matching::perfect() const { inline bool Graph_matching::multi_augment() { if (perfect()) return false; - Layered_neighbors_finder layered_nf = *layering(); + Layered_neighbors_finder layered_nf(layering()); int max_depth = layered_nf.vlayers_number()*2 - 1; double rn = sqrt(G::size()); // verification of a necessary criterion in order to shortcut if possible @@ -103,7 +101,7 @@ inline void Graph_matching::set_r(double r) { inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) { //V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices. - std::deque path; + std::vector path; path.emplace_back(u_start_index); do { if (static_cast(path.size()) > max_depth) { @@ -129,19 +127,19 @@ inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u return true; } -inline std::shared_ptr Graph_matching::layering() const { +inline Layered_neighbors_finder Graph_matching::layering() const { std::list u_vertices(unmatched_in_u); std::list v_vertices; Neighbors_finder nf(r); for (int v_point_index = 0; v_point_index < G::size(); ++v_point_index) nf.add(v_point_index); - std::shared_ptr layered_nf(new Layered_neighbors_finder(r)); + Layered_neighbors_finder layered_nf(r); for(int layer = 0; !u_vertices.empty(); layer++) { // one layer is one step in the BFS for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) { - std::shared_ptr> u_succ(nf.pull_all_near(*it1)); - for (auto it2 = u_succ->begin(); it2 != u_succ->end(); ++it2) { - layered_nf->add(*it2, layer); + std::vector u_succ(nf.pull_all_near(*it1)); + for (auto it2 = u_succ.begin(); it2 != u_succ.end(); ++it2) { + layered_nf.add(*it2, layer); v_vertices.emplace_back(*it2); } } @@ -162,7 +160,7 @@ inline std::shared_ptr Graph_matching::layering() cons return layered_nf; } -inline void Graph_matching::update(std::deque& path) { +inline void Graph_matching::update(std::vector& path) { unmatched_in_u.remove(path.front()); for (auto it = path.cbegin(); it != path.cend(); ++it) { // Be careful, the iterator is incremented twice each time diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h index 1bd5879c..e8560559 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -46,7 +46,7 @@ public: /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ int pull_near(int u_point_index); /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */ - std::shared_ptr< std::list > pull_all_near(int u_point_index); + std::vector pull_all_near(int u_point_index); private: const double r; @@ -118,11 +118,11 @@ inline int Neighbors_finder::pull_near(int u_point_index) { return tmp; } -inline std::shared_ptr< std::list > Neighbors_finder::pull_all_near(int u_point_index) { - std::shared_ptr< std::list > all_pull(planar_neighbors_f.pull_all_near(u_point_index)); +inline std::vector Neighbors_finder::pull_all_near(int u_point_index) { + std::vector all_pull(planar_neighbors_f.pull_all_near(u_point_index)); int last_pull = pull_near(u_point_index); while (last_pull != null_point_index()) { - all_pull->emplace_back(last_pull); + all_pull.emplace_back(last_pull); last_pull = pull_near(u_point_index); } return all_pull; diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h index 1f6d6683..7af149e4 100644 --- a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h @@ -25,11 +25,6 @@ #include #include -#include -#include -#include -#include -#include #include namespace Gudhi { @@ -60,7 +55,7 @@ public: /** \internal \brief Returns size = |U| = |V|. */ static int size(); /** \internal \brief Returns the O(n^2) sorted distances between the points. */ - static std::shared_ptr< std::vector > sorted_distances(); + static std::vector sorted_distances(); /** \internal \brief Returns an upper bound of the diameter of the convex hull */ static double diameter(); @@ -126,15 +121,14 @@ inline int G::size() { return static_cast (u.size() + v.size()); } -inline std::shared_ptr< std::vector > G::sorted_distances() { +inline std::vector G::sorted_distances() { // could be optimized std::set sorted_distances; sorted_distances.emplace(0.); for (int u_point_index = 0; u_point_index < size(); ++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)); - std::shared_ptr< std::vector > sd_up(new std::vector(sorted_distances.cbegin(), sorted_distances.cend())); - return sd_up; + return std::vector(sorted_distances.begin(),sorted_distances.end()); } inline Internal_point G::get_u_point(int u_point_index) { diff --git a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h index 8bd21267..f574231e 100644 --- a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h @@ -23,9 +23,6 @@ #ifndef PLANAR_NEIGHBORS_FINDER_H_ #define PLANAR_NEIGHBORS_FINDER_H_ -#include -#include - // Inclusion order is important for CGAL patch #include "CGAL/Kd_tree_node.h" #include "CGAL/Kd_tree.h" @@ -62,7 +59,7 @@ public: /** \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. */ int pull_near(int u_point_index); /** \internal \brief Provide and remove all the V points near to the U point given as parameter. */ - virtual std::shared_ptr< std::list > pull_all_near(int u_point_index); + std::vector pull_all_near(int u_point_index); private: double r; @@ -91,7 +88,7 @@ public: /** \internal \brief Provide a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ int pull_near(int u_point_index); /** \internal \brief Provide and remove all the V points near to the U point given as parameter. */ - virtual std::shared_ptr< std::list > pull_all_near(int u_point_index); + virtual std::vector pull_all_near(int u_point_index); private: double r; @@ -148,8 +145,8 @@ inline int Naive_pnf::pull_near(int u_point_index) { return null_point_index(); } -inline std::shared_ptr< std::list > Naive_pnf::pull_all_near(int u_point_index) { - std::shared_ptr< std::list > all_pull(new std::list); +inline std::vector Naive_pnf::pull_all_near(int u_point_index) { + std::vector all_pull; Internal_point u_point = G::get_u_point(u_point_index); int i0 = static_cast(u_point.x()/r); int j0 = static_cast(u_point.y()/r); @@ -157,7 +154,7 @@ inline std::shared_ptr< std::list > Naive_pnf::pull_all_near(int u_point_in for(int j = 1; j<= 3; j++) for(auto it = grid.find(std::make_pair(i0 +(i%3)-1, j0+(j%3)-1)); it!=grid.end();) if (G::distance(u_point_index, it->second) <= r) { - all_pull->emplace_back(it->second); + all_pull.emplace_back(it->second); it = grid.erase(it); } else it++; @@ -206,11 +203,11 @@ inline int Cgal_pnf::pull_near(int u_point_index){ return tmp; } -inline std::shared_ptr< std::list > Cgal_pnf::pull_all_near(int u_point_index) { - std::shared_ptr< std::list > all_pull(new std::list); +inline std::vector Cgal_pnf::pull_all_near(int u_point_index) { + std::vector all_pull; int last_pull = pull_near(u_point_index); while (last_pull != null_point_index()) { - all_pull->emplace_back(last_pull); + all_pull.emplace_back(last_pull); last_pull = pull_near(u_point_index); } return all_pull; -- cgit v1.2.3 From 5bd04a9433462965aecf000e5044e70405e968ff Mon Sep 17 00:00:00 2001 From: fgodi Date: Wed, 23 Nov 2016 16:57:49 +0000 Subject: fix git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1780 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 9001cd1b33dfea3b34e5d12ef8e20b5624a50ba9 --- .../concept/Persistence_diagram.h | 11 +++++------ .../example/bottleneck_example.cpp | 2 +- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 20 +++++++------------- .../include/gudhi/Persistence_diagrams_graph.h | 4 ++-- .../test/bottleneck_unit_test.cpp | 1 + 5 files changed, 16 insertions(+), 22 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/concept/Persistence_diagram.h b/src/Bottleneck_distance/concept/Persistence_diagram.h index 4c69343c..1ab768ba 100644 --- a/src/Bottleneck_distance/concept/Persistence_diagram.h +++ b/src/Bottleneck_distance/concept/Persistence_diagram.h @@ -27,23 +27,22 @@ namespace Gudhi { namespace bottleneck_distance { -/** \brief Concept of persistence diagram point. The double first is the birth of the component and the double second is the death of the component. +/** \brief Concept of persistence diagram point. get<0>() must return the birth of the component and get<1>() its death. * * \ingroup bottleneck_distance */ struct Diagram_point{ - double first; - double second; + double get(); }; -/** \brief Concept of persistence diagram. +/** \brief Concept of persistence diagram. It's a range of Diagram_point. * * \ingroup bottleneck_distance */ struct Persistence_Diagram { - const_iterator cbegin() const; - const_iterator cend() const; + const_iterator begin(); + const_iterator end(); }; } // namespace bottleneck_distance diff --git a/src/Bottleneck_distance/example/bottleneck_example.cpp b/src/Bottleneck_distance/example/bottleneck_example.cpp index fd0f8ed5..b1b98a82 100644 --- a/src/Bottleneck_distance/example/bottleneck_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_example.cpp @@ -73,5 +73,5 @@ int main( int argc , char** argv ) tolerance = atof( argv[3] ); } double b = Gudhi::bottleneck_distance::compute(diag1, diag2, tolerance); - std::cout << "The distance between the diagrams is : " << b << ". The tolerace is : " << tolerance << std::endl; + std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; } diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index 1ae7788c..71845e25 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -29,17 +29,6 @@ namespace Gudhi { namespace bottleneck_distance { -/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams. You get an additive e-approximation. - * - * - * \ingroup bottleneck_distance - */ -template -double compute(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e = 0.); - -template -double compute_exactly(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2); - template double compute_exactly(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2) { G::initialize(diag1, diag2, 0.); @@ -66,9 +55,14 @@ double compute_exactly(const Persistence_diagram1 &diag1, const Persistence_diag return sd.at(idmin); } +/** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams. + * If the last parameter e is not 0 (default value if not explicited), you get an additive e-approximation. + * + * \ingroup bottleneck_distance + */ template -double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { - if(e< std::numeric_limits::min()) +double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=0.) { + if(e == 0.) return compute_exactly(diag1, diag2); G::initialize(diag1, diag2, e); int in = G::diameter()/e; diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h index 7af149e4..7d37b9e0 100644 --- a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h @@ -83,10 +83,10 @@ inline void G::initialize(const Persistence_diagram1 &diag1, v.clear(); for (auto it = diag1.cbegin(); it != diag1.cend(); ++it) if (it->second - it->first > e) - u.push_back(Internal_point(it->first, it->second, u.size())); + 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 - it->first > e) - v.push_back(Internal_point(it->first, it->second, v.size())); + v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size())); if (u.size() < v.size()) swap(u, v); } diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp index b757ce54..e2cd3c05 100644 --- a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp +++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp @@ -186,4 +186,5 @@ 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.); } -- cgit v1.2.3 From e15d7f55c337596e988cd84425bbfc4815913074 Mon Sep 17 00:00:00 2001 From: fgodi Date: Fri, 25 Nov 2016 13:17:50 +0000 Subject: copyrigths git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1783 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: c045f59d4f2810aa2193add84ac97a2632888dc4 --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 6 +- .../include/gudhi/Construct_coord_iterator.h | 4 +- .../include/gudhi/Graph_matching.h | 4 +- .../include/gudhi/Internal_point.h | 4 +- .../include/gudhi/Neighbors_finder.h | 4 +- .../include/gudhi/Persistence_diagrams_graph.h | 7 +- .../include/gudhi/Planar_neighbors_finder.h | 89 +++------------------- 7 files changed, 26 insertions(+), 92 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index 71845e25..bb0a13d2 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Francois Godi + * Author: Francois Godi * - * Copyright (C) 2015 INRIA Sophia-Antipolis (France) + * Copyright (C) 2015 INRIA (France) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -65,7 +65,7 @@ double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &di if(e == 0.) return compute_exactly(diag1, diag2); G::initialize(diag1, diag2, e); - int in = G::diameter()/e; + int in = G::diameter()/e + 1; int idmin = 0; int idmax = in; // alpha can be modified, this will change the complexity diff --git a/src/Bottleneck_distance/include/gudhi/Construct_coord_iterator.h b/src/Bottleneck_distance/include/gudhi/Construct_coord_iterator.h index 278b7955..0959dd03 100644 --- a/src/Bottleneck_distance/include/gudhi/Construct_coord_iterator.h +++ b/src/Bottleneck_distance/include/gudhi/Construct_coord_iterator.h @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Francois Godi + * Author: Francois Godi * - * Copyright (C) 2015 INRIA Sophia-Antipolis (France) + * Copyright (C) 2015 INRIA (France) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h index c50c3a9e..fa20b2a2 100644 --- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Francois Godi + * Author: Francois Godi * - * Copyright (C) 2015 INRIA Sophia-Antipolis (France) + * Copyright (C) 2015 INRIA (France) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h index 90b5dd14..78aad470 100644 --- a/src/Bottleneck_distance/include/gudhi/Internal_point.h +++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Francois Godi + * Author: Francois Godi * - * Copyright (C) 2015 INRIA Sophia-Antipolis (France) + * Copyright (C) 2015 INRIA (France) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h index e8560559..20b47d0b 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Francois Godi + * Author: Francois Godi * - * Copyright (C) 2015 INRIA Sophia-Antipolis (France) + * Copyright (C) 2015 INRIA (France) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h index 7d37b9e0..8242ce2b 100644 --- a/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h +++ b/src/Bottleneck_distance/include/gudhi/Persistence_diagrams_graph.h @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Francois Godi + * Author: Francois Godi * - * Copyright (C) 2015 INRIA Sophia-Antipolis (France) + * Copyright (C) 2015 INRIA (France) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -65,8 +65,7 @@ private: static Internal_point get_u_point(int u_point_index); static Internal_point get_v_point(int v_point_index); - friend class Naive_pnf; - friend class Cgal_pnf; + friend class Planar_neighbors_finder; }; /** \internal \typedef \brief Shorter alias */ diff --git a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h index f574231e..c1f4d908 100644 --- a/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Planar_neighbors_finder.h @@ -2,9 +2,9 @@ * (Geometric Understanding in Higher Dimensions) is a generic C++ * library for computational topology. * - * Author(s): Francois Godi + * Author: Francois Godi * - * Copyright (C) 2015 INRIA Sophia-Antipolis (France) + * Copyright (C) 2015 INRIA (France) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -67,7 +67,7 @@ private: std::multimap,int> grid; }; -class Cgal_pnf { +class Planar_neighbors_finder { typedef CGAL::Dimension_tag<2> D; typedef CGAL::Search_traits Traits; @@ -78,7 +78,7 @@ class Cgal_pnf { public: /** \internal \brief Constructor taking the near distance definition as parameter. */ - Cgal_pnf(double r_); + Planar_neighbors_finder(double r_); /** \internal \brief A point added will be possibly pulled. */ void add(int v_point_index); /** \internal \brief A point manually removed will no longer be possibly pulled. */ @@ -96,79 +96,14 @@ private: Kd_tree kd_t; }; -/** \internal \typedef \brief Planar_neighbors_finder is the used implementation. */ -typedef Cgal_pnf Planar_neighbors_finder; - -inline Naive_pnf::Naive_pnf(double r_) : - r(r_), grid() { } - - -inline std::pair Naive_pnf::get_v_key(int v_point_index) const{ - Internal_point v_point = G::get_v_point(v_point_index); - return std::make_pair(static_cast(v_point.x()/r), static_cast(v_point.y()/r)); -} - -inline void Naive_pnf::add(int v_point_index) { - grid.emplace(get_v_key(v_point_index),v_point_index); -} - -inline void Naive_pnf::remove(int v_point_index) { - if(v_point_index != null_point_index()) - for(auto it = grid.find(get_v_key(v_point_index)); it!=grid.end(); it++) - if(it->second==v_point_index){ - grid.erase(it); - return; - } -} - -inline bool Naive_pnf::contains(int v_point_index) const { - if(v_point_index == null_point_index()) - return false; - for(auto it = grid.find(get_v_key(v_point_index)); it!=grid.end(); it++) - if(it->second==v_point_index) - return true; - return false; -} - -inline int Naive_pnf::pull_near(int u_point_index) { - Internal_point u_point = G::get_u_point(u_point_index); - int i0 = static_cast(u_point.x()/r); - int j0 = static_cast(u_point.y()/r); - for(int i = 1; i<= 3; i++) - for(int j = 1; j<= 3; j++) - for(auto it = grid.find(std::make_pair(i0 +(i%3)-1, j0+(j%3)-1)); it!=grid.end(); it++) - if (G::distance(u_point_index, it->second) <= r) { - int tmp = it->second; - grid.erase(it); - return tmp; - } - return null_point_index(); -} - -inline std::vector Naive_pnf::pull_all_near(int u_point_index) { - std::vector all_pull; - Internal_point u_point = G::get_u_point(u_point_index); - int i0 = static_cast(u_point.x()/r); - int j0 = static_cast(u_point.y()/r); - for(int i = 1; i<= 3; i++) - for(int j = 1; j<= 3; j++) - for(auto it = grid.find(std::make_pair(i0 +(i%3)-1, j0+(j%3)-1)); it!=grid.end();) - if (G::distance(u_point_index, it->second) <= r) { - all_pull.emplace_back(it->second); - it = grid.erase(it); - } - else it++; - return all_pull; -} - /** \internal \brief Constructor taking the near distance definition as parameter. */ -inline Cgal_pnf::Cgal_pnf(double r_) +inline Planar_neighbors_finder::Planar_neighbors_finder(double r_) : r(r_), contents(), kd_t() {} /** \internal \brief A point added will be possibly pulled. */ -inline void Cgal_pnf::add(int v_point_index){ +inline void Planar_neighbors_finder::add(int v_point_index){ if(v_point_index == null_point_index()) return; contents.insert(v_point_index); @@ -176,22 +111,20 @@ inline void Cgal_pnf::add(int v_point_index){ } /** \internal \brief A point manually removed will no longer be possibly pulled. */ -inline void Cgal_pnf::remove(int v_point_index){ - if(contains(v_point_index)){ +inline void Planar_neighbors_finder::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 ? */ -inline bool Cgal_pnf::contains(int v_point_index) const{ +inline bool Planar_neighbors_finder::contains(int v_point_index) const{ if(v_point_index == null_point_index()) return false; return contents.count(v_point_index)>0; } /** \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){ +inline int Planar_neighbors_finder::pull_near(int u_point_index){ Internal_point u_point = G::get_u_point(u_point_index); std::vector w = {1., 1.}; K_neighbor_search search(kd_t, u_point, 0., true, Distance(0, 2, w)); @@ -199,11 +132,13 @@ inline int Cgal_pnf::pull_near(int u_point_index){ if(it==search.end() || G::distance(u_point_index, it->first.point_index) > r) return null_point_index(); int tmp = it->first.point_index; + if(!contains(tmp)) + std::cout << "!! A kd_tree returns a point (Point_index:" << tmp << ") previously removed !!" << std::endl; remove(tmp); return tmp; } -inline std::vector Cgal_pnf::pull_all_near(int u_point_index) { +inline std::vector Planar_neighbors_finder::pull_all_near(int u_point_index) { std::vector all_pull; int last_pull = pull_near(u_point_index); while (last_pull != null_point_index()) { -- cgit v1.2.3 From cc1e554f26765bf9993079bef608d9bc4a3308c8 Mon Sep 17 00:00:00 2001 From: fgodi Date: Mon, 28 Nov 2016 14:15:22 +0000 Subject: graph no longer static git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1790 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 6bbb8dfaffd09a81ccafac6919e073ed74b9c7e6 --- .../concept/Persistence_diagram.h | 4 +- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 22 +++--- .../include/gudhi/Construct_coord_iterator.h | 42 ---------- .../include/gudhi/Graph_matching.h | 18 +++-- .../include/gudhi/Internal_point.h | 18 ++++- .../include/gudhi/Neighbors_finder.h | 36 +++++---- .../test/bottleneck_unit_test.cpp | 92 ++++++++++++---------- 7 files changed, 107 insertions(+), 125 deletions(-) delete mode 100644 src/Bottleneck_distance/include/gudhi/Construct_coord_iterator.h (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/concept/Persistence_diagram.h b/src/Bottleneck_distance/concept/Persistence_diagram.h index fe13e859..2ca10094 100644 --- a/src/Bottleneck_distance/concept/Persistence_diagram.h +++ b/src/Bottleneck_distance/concept/Persistence_diagram.h @@ -41,8 +41,8 @@ struct Diagram_point{ */ struct Persistence_Diagram { - const_iterator begin(); - const_iterator end(); + iterator begin(); + iterator end(); }; } // namespace bottleneck_distance diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index bb0a13d2..52de30be 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -31,14 +31,14 @@ namespace bottleneck_distance { template double compute_exactly(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2) { - G::initialize(diag1, diag2, 0.); - std::vector sd(G::sorted_distances()); + Persistence_graph g(diag1, diag2, 0.); + std::vector sd(g.sorted_distances()); int idmin = 0; int idmax = sd.size() - 1; - // alpha can be modified, this will change the complexity + // alpha can change the complexity double alpha = pow(sd.size(), 0.25); - Graph_matching m; - Graph_matching biggest_unperfect; + Graph_matching m(g); + Graph_matching biggest_unperfect(g); while (idmin != idmax) { int step = static_cast((idmax - idmin) / alpha); m.set_r(sd.at(idmin + step)); @@ -64,14 +64,14 @@ template double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=0.) { if(e == 0.) return compute_exactly(diag1, diag2); - G::initialize(diag1, diag2, e); - int in = G::diameter()/e + 1; + Persistence_graph g(diag1, diag2, e); + int in = g.diameter_bound()/e + 1; int idmin = 0; int idmax = in; - // alpha can be modified, this will change the complexity - double alpha = pow(in, 0.25); - Graph_matching m; - Graph_matching biggest_unperfect; + // alpha can change the complexity + double alpha = pow(in, 0.10); + Graph_matching m(g); + Graph_matching biggest_unperfect(g); while (idmin != idmax) { int step = static_cast((idmax - idmin) / alpha); m.set_r(e*(idmin + step)); diff --git a/src/Bottleneck_distance/include/gudhi/Construct_coord_iterator.h b/src/Bottleneck_distance/include/gudhi/Construct_coord_iterator.h deleted file mode 100644 index 0959dd03..00000000 --- a/src/Bottleneck_distance/include/gudhi/Construct_coord_iterator.h +++ /dev/null @@ -1,42 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author: Francois Godi - * - * Copyright (C) 2015 INRIA (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#ifndef SRC_BOTTLENECK_INCLUDE_CGAL_MISCELLANEOUS_H_ -#define SRC_BOTTLENECK_INCLUDE_CGAL_MISCELLANEOUS_H_ - -#include - -namespace CGAL { - -typedef Gudhi::bottleneck_distance::Internal_point Internal_point; - -struct Construct_coord_iterator { - typedef const double* result_type; - const double* operator()(const Internal_point& p) const - { return p.vec; } - const double* operator()(const Internal_point& p, int) const - { return p.vec+2; } -}; - -} //namespace CGAL - -#endif // SRC_BOTTLENECK_INCLUDE_CGAL_MISCELLANEOUS_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h index fa20b2a2..10b7073a 100644 --- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -36,7 +36,7 @@ namespace bottleneck_distance { class Graph_matching { public: /** \internal \brief Constructor constructing an empty matching. */ - explicit Graph_matching(); + explicit Graph_matching(Persistence_graph &g); /** \internal \brief Copy operator. */ Graph_matching& operator=(const Graph_matching& m); /** \internal \brief Is the matching perfect ? */ @@ -47,6 +47,7 @@ public: void set_r(double r); private: + Persistence_graph& g; double r; /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */ std::vector v_to_u; @@ -61,13 +62,14 @@ private: void update(std::vector & path); }; -inline Graph_matching::Graph_matching() - : r(0.), v_to_u(G::size(), null_point_index()), unmatched_in_u() { - for (int u_point_index = 0; u_point_index < G::size(); ++u_point_index) +inline Graph_matching::Graph_matching(Persistence_graph& g) + : g(g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u() { + for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index) unmatched_in_u.emplace_back(u_point_index); } inline Graph_matching& Graph_matching::operator=(const Graph_matching& m) { + g = m.g; r = m.r; v_to_u = m.v_to_u; unmatched_in_u = m.unmatched_in_u; @@ -83,7 +85,7 @@ inline bool Graph_matching::multi_augment() { return false; Layered_neighbors_finder layered_nf(layering()); int max_depth = layered_nf.vlayers_number()*2 - 1; - double rn = sqrt(G::size()); + double rn = sqrt(g.size()); // verification of a necessary criterion in order to shortcut if possible if (max_depth <0 || (unmatched_in_u.size() > rn && max_depth >= rn)) return false; @@ -130,10 +132,10 @@ inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u inline Layered_neighbors_finder Graph_matching::layering() const { std::list u_vertices(unmatched_in_u); std::list v_vertices; - Neighbors_finder nf(r); - for (int v_point_index = 0; v_point_index < G::size(); ++v_point_index) + Neighbors_finder nf(g, r); + for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index) nf.add(v_point_index); - Layered_neighbors_finder layered_nf(r); + Layered_neighbors_finder layered_nf(g, r); for(int layer = 0; !u_vertices.empty(); layer++) { // one layer is one step in the BFS for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) { diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h index 78aad470..1a050d48 100644 --- a/src/Bottleneck_distance/include/gudhi/Internal_point.h +++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h @@ -35,7 +35,7 @@ struct Internal_point { double vec[2]; int 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; } + Internal_point(double x, double y, int p_i) { 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 ]; } @@ -44,7 +44,7 @@ struct Internal_point { { return point_index==p.point_index; } - bool operator!=(const Internal_point& p) const { return ! (*this == p); } + bool operator!=(const Internal_point& p) const { return !(*this == p); } }; inline int null_point_index() { @@ -55,4 +55,18 @@ inline int null_point_index() { } // namespace Gudhi +namespace CGAL { + +typedef Gudhi::bottleneck_distance::Internal_point Internal_point; + +struct Construct_coord_iterator { + typedef const double* result_type; + const double* operator()(const Internal_point& p) const + { return p.vec; } + const double* operator()(const Internal_point& p, int) const + { return p.vec+2; } +}; + +} //namespace CGAL + #endif // INTERNAL_POINT_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h index f2d9abb2..f28e21a2 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -31,8 +31,8 @@ #include #include -#include -#include +#include +#include #include @@ -57,7 +57,7 @@ class Neighbors_finder { public: /** \internal \brief Constructor taking the near distance definition as parameter. */ - Neighbors_finder(double r); + Neighbors_finder(const Persistence_graph& g, double r); /** \internal \brief A point added will be possibly pulled. */ void add(int v_point_index); /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ @@ -66,6 +66,7 @@ public: std::vector pull_all_near(int u_point_index); private: + const Persistence_graph& g; const double r; Kd_tree kd_t; std::unordered_set projections_f; @@ -81,7 +82,7 @@ private: class Layered_neighbors_finder { public: /** \internal \brief Constructor taking the near distance definition as parameter. */ - Layered_neighbors_finder(double r); + Layered_neighbors_finder(const Persistence_graph& g, double r); /** \internal \brief A point added will be possibly pulled. */ void add(int v_point_index, int vlayer); /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ @@ -90,43 +91,44 @@ public: int vlayers_number() const; private: + const Persistence_graph& g; const double r; std::vector> neighbors_finder; }; -inline Neighbors_finder::Neighbors_finder(double r) : - r(r), kd_t(), projections_f() { } +inline Neighbors_finder::Neighbors_finder(const Persistence_graph& g, double r) : + g(g), r(r), kd_t(), projections_f() { } inline void Neighbors_finder::add(int v_point_index) { - if (G::on_the_v_diagonal(v_point_index)) + if (g.on_the_v_diagonal(v_point_index)) projections_f.emplace(v_point_index); else - kd_t.insert(G::get_v_point(v_point_index)); + kd_t.insert(g.get_v_point(v_point_index)); } inline int Neighbors_finder::pull_near(int u_point_index) { int tmp; - int c = G::corresponding_point_in_v(u_point_index); - if (G::on_the_u_diagonal(u_point_index) && !projections_f.empty()){ + int c = g.corresponding_point_in_v(u_point_index); + if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()){ //Any pair of projection is at distance 0 tmp = *projections_f.cbegin(); projections_f.erase(tmp); } - else if (projections_f.count(c) && (G::distance(u_point_index, c) <= r)){ + else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)){ //Is the query point near to its projection ? tmp = c; projections_f.erase(tmp); } else{ //Is the query point near to a V point in the plane ? - Internal_point u_point = G::get_u_point(u_point_index); + Internal_point u_point = g.get_u_point(u_point_index); std::vector 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) + if(it==search.end() || g.distance(u_point_index, it->first.point_index) > r) return null_point_index(); tmp = it->first.point_index; - kd_t.remove(G::get_v_point(tmp)); + kd_t.remove(g.get_v_point(tmp)); } return tmp; } @@ -141,12 +143,12 @@ inline std::vector Neighbors_finder::pull_all_near(int u_point_index) { return all_pull; } -inline Layered_neighbors_finder::Layered_neighbors_finder(double r) : - r(r), neighbors_finder() { } +inline Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_graph& g, double r) : + g(g), r(r), neighbors_finder() { } inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) { for (int l = neighbors_finder.size(); l <= vlayer; l++) - neighbors_finder.emplace_back(std::shared_ptr(new Neighbors_finder(r))); + neighbors_finder.emplace_back(std::shared_ptr(new Neighbors_finder(g, r))); neighbors_finder.at(vlayer)->add(v_point_index); } diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp index 7edfa062..2237b073 100644 --- a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp +++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp @@ -34,11 +34,13 @@ int n1 = 81; // a natural number >0 int n2 = 180; // a natural number >0 double upper_bound = 406.43; // any real >0 + +std::uniform_real_distribution unif(0.,upper_bound); +std::default_random_engine re; +std::vector< std::pair > v1, v2; + BOOST_AUTO_TEST_CASE(persistence_diagrams_graph){ // Random construction - std::uniform_real_distribution unif(0.,upper_bound); - std::default_random_engine re; - std::vector< std::pair > v1, v2; for (int i = 0; i < n1; i++) { double a = unif(re); double b = unif(re); @@ -49,83 +51,87 @@ BOOST_AUTO_TEST_CASE(persistence_diagrams_graph){ double b = unif(re); v2.emplace_back(std::min(a,b), std::max(a,b)); } - G::initialize(v1, v2, 0.); - std::vector d(G::sorted_distances()); + Persistence_graph g(v1, v2, 0.); + std::vector d(g.sorted_distances()); // - BOOST_CHECK(!G::on_the_u_diagonal(n1-1)); - BOOST_CHECK(!G::on_the_u_diagonal(n1)); - BOOST_CHECK(!G::on_the_u_diagonal(n2-1)); - BOOST_CHECK(G::on_the_u_diagonal(n2)); - BOOST_CHECK(!G::on_the_v_diagonal(n1-1)); - BOOST_CHECK(G::on_the_v_diagonal(n1)); - BOOST_CHECK(G::on_the_v_diagonal(n2-1)); - BOOST_CHECK(G::on_the_v_diagonal(n2)); + BOOST_CHECK(!g.on_the_u_diagonal(n1-1)); + BOOST_CHECK(!g.on_the_u_diagonal(n1)); + BOOST_CHECK(!g.on_the_u_diagonal(n2-1)); + BOOST_CHECK(g.on_the_u_diagonal(n2)); + BOOST_CHECK(!g.on_the_v_diagonal(n1-1)); + BOOST_CHECK(g.on_the_v_diagonal(n1)); + BOOST_CHECK(g.on_the_v_diagonal(n2-1)); + BOOST_CHECK(g.on_the_v_diagonal(n2)); // - BOOST_CHECK(G::corresponding_point_in_u(0)==n2); - BOOST_CHECK(G::corresponding_point_in_u(n1)==0); - BOOST_CHECK(G::corresponding_point_in_v(0)==n1); - BOOST_CHECK(G::corresponding_point_in_v(n2)==0); + BOOST_CHECK(g.corresponding_point_in_u(0)==n2); + BOOST_CHECK(g.corresponding_point_in_u(n1)==0); + BOOST_CHECK(g.corresponding_point_in_v(0)==n1); + BOOST_CHECK(g.corresponding_point_in_v(n2)==0); // - BOOST_CHECK(G::size()==(n1+n2)); + BOOST_CHECK(g.size()==(n1+n2)); // BOOST_CHECK((int) d.size() <= (n1+n2)*(n1+n2) - n1*n2 + 1); - 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(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_AUTO_TEST_CASE(neighbors_finder) { - Neighbors_finder nf(1.); + Persistence_graph g(v1, v2, 0.); + Neighbors_finder nf(g, 1.); for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2) nf.add(v_point_index); // int v_point_index_1 = nf.pull_near(n2/2); - BOOST_CHECK((v_point_index_1 == -1) || (G::distance(n2/2,v_point_index_1)<=1.)); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); std::vector l = nf.pull_all_near(n2/2); bool v = true; for(auto it = l.cbegin(); it != l.cend(); ++it) - v = v && (G::distance(n2/2,*it)>1.); + v = v && (g.distance(n2/2,*it)>1.); BOOST_CHECK(v); int v_point_index_2 = nf.pull_near(n2/2); BOOST_CHECK(v_point_index_2 == -1); } BOOST_AUTO_TEST_CASE(layered_neighbors_finder) { - Layered_neighbors_finder lnf(1.); + Persistence_graph g(v1, v2, 0.); + Layered_neighbors_finder lnf(g, 1.); for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2) lnf.add(v_point_index, v_point_index % 7); // int v_point_index_1 = lnf.pull_near(n2/2,6); - BOOST_CHECK((v_point_index_1 == -1) || (G::distance(n2/2,v_point_index_1)<=1.)); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); int v_point_index_2 = lnf.pull_near(n2/2,6); BOOST_CHECK(v_point_index_2 == -1); v_point_index_1 = lnf.pull_near(n2/2,0); - BOOST_CHECK((v_point_index_1 == -1) || (G::distance(n2/2,v_point_index_1)<=1.)); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); v_point_index_2 = lnf.pull_near(n2/2,0); BOOST_CHECK(v_point_index_2 == -1); } BOOST_AUTO_TEST_CASE(graph_matching) { - Graph_matching m1; + Persistence_graph g(v1, v2, 0.); + Graph_matching m1(g); m1.set_r(0.); int e = 0; while (m1.multi_augment()) ++e; + BOOST_CHECK(e > 0); BOOST_CHECK(e <= 2*sqrt(2*(n1+n2))); Graph_matching m2 = m1; BOOST_CHECK(!m2.multi_augment()); -- cgit v1.2.3 From 7449716dca77dd81759d024ae1e9dbfcfeb202e7 Mon Sep 17 00:00:00 2001 From: fgodi Date: Mon, 28 Nov 2016 15:25:17 +0000 Subject: compute and compute_exactly merged git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1792 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 7c4e66f5449759f7d99319ca14989b2d74e96e20 --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 62 +++------------------- .../include/gudhi/Neighbors_finder.h | 4 +- 2 files changed, 9 insertions(+), 57 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index 52de30be..59df655e 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -29,32 +29,6 @@ namespace Gudhi { namespace bottleneck_distance { -template -double compute_exactly(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2) { - Persistence_graph g(diag1, diag2, 0.); - std::vector sd(g.sorted_distances()); - int idmin = 0; - int idmax = sd.size() - 1; - // alpha can change the complexity - double alpha = pow(sd.size(), 0.25); - Graph_matching m(g); - Graph_matching biggest_unperfect(g); - while (idmin != idmax) { - int step = static_cast((idmax - idmin) / alpha); - m.set_r(sd.at(idmin + step)); - while (m.multi_augment()); - //The above while compute a maximum matching (according to the r setted before) - if (m.perfect()) { - idmax = idmin + step; - m = biggest_unperfect; - } else { - biggest_unperfect = m; - idmin = idmin + step + 1; - } - } - return sd.at(idmin); -} - /** \brief Function to use in order to compute the Bottleneck distance between two persistence diagrams. * If the last parameter e is not 0 (default value if not explicited), you get an additive e-approximation. * @@ -62,19 +36,19 @@ double compute_exactly(const Persistence_diagram1 &diag1, const Persistence_diag */ template double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=0.) { - if(e == 0.) - return compute_exactly(diag1, diag2); Persistence_graph g(diag1, diag2, e); - int in = g.diameter_bound()/e + 1; + std::vector sd; + if(e == 0.) + sd = g.sorted_distances(); int idmin = 0; - int idmax = in; + int idmax = e==0. ? sd.size() - 1 : g.diameter_bound()/e + 1; // alpha can change the complexity - double alpha = pow(in, 0.10); + 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); - m.set_r(e*(idmin + step)); + 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) if (m.perfect()) { @@ -85,7 +59,7 @@ double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &di idmin = idmin + step + 1; } } - return e*(idmin); + return e == 0. ? sd.at(idmin) : e*(idmin); } @@ -95,25 +69,3 @@ double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &di } // namespace Gudhi #endif // BOTTLENECK_H_ - -/* Dichotomic version -template -double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) { - if(e< std::numeric_limits::min()) - return compute_exactly(diag1, diag2); - G::initialize(diag1, diag2, e); - double d = 0.; - double f = G::diameter(); - while (f-d > e){ - Graph_matching m; - m.set_r((d+f)/2.); - while (m.multi_augment()); - //The above while compute a maximum matching (according to the r setted before) - if (m.perfect()) - f = (d+f)/2.; - else - d= (d+f)/2.; - } - return d; -} */ - diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h index f28e21a2..90ddabef 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -122,8 +122,8 @@ inline int Neighbors_finder::pull_near(int u_point_index) { else{ //Is the query point near to a V point in the plane ? Internal_point u_point = g.get_u_point(u_point_index); - std::vector w = {1., 1.}; - K_neighbor_search search(kd_t, u_point, 0., true, Distance(0, 2, w)); + std::array w = { {1., 1.} }; + K_neighbor_search search(kd_t, u_point, 0., true, Distance(0, 2, w.begin(), w.end())); auto it = search.begin(); if(it==search.end() || g.distance(u_point_index, it->first.point_index) > r) return null_point_index(); -- cgit v1.2.3 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(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') 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 From d0db67f8ff185a7698d2f6643d86e43955bab882 Mon Sep 17 00:00:00 2001 From: fgodi Date: Tue, 6 Dec 2016 16:20:10 +0000 Subject: 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 --- .../concept/Persistence_diagram.h | 16 ++++------- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 11 +++++--- .../include/gudhi/Persistence_graph.h | 33 +++++++++++++--------- 3 files changed, 33 insertions(+), 27 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') 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::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. + * Birth must be 0 or positive, death should be larger than birth, death can be std::numeric_limits::infinity() for components which stay alive. * * \ingroup bottleneck_distance */ -struct Diagram_point{ - double get(); -}; +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 begin(); - iterator 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 +#include 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 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::infinity()) return std::numeric_limits::infinity(); std::vector 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 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 u; std::vector v; - int alive_count; + std::vector u_alive; + std::vector v_alive; }; template 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::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::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 (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::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 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::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::infinity() ? it->x() : it->y()); + max = std::max(max, it->y()); return max; } + } // namespace bottleneck_distance } // namespace Gudhi -- cgit v1.2.3 From 9906477eb45f16cfafe7186482f4c44606c575dc Mon Sep 17 00:00:00 2001 From: fgodi Date: Thu, 15 Dec 2016 10:00:19 +0000 Subject: namespace git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1877 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 0d15f01525272ec8bea7436fa18dbac213becf0a --- src/Bottleneck_distance/example/bottleneck_basic_example.cpp | 4 ++-- src/Bottleneck_distance/example/bottleneck_read_file_example.cpp | 2 +- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 6 +++--- src/Bottleneck_distance/include/gudhi/Graph_matching.h | 4 ++-- src/Bottleneck_distance/include/gudhi/Internal_point.h | 4 ++-- src/Bottleneck_distance/include/gudhi/Neighbors_finder.h | 4 ++-- src/Bottleneck_distance/include/gudhi/Persistence_graph.h | 4 ++-- src/Bottleneck_distance/test/bottleneck_chrono.cpp | 6 +++--- src/Bottleneck_distance/test/bottleneck_unit_test.cpp | 8 ++++---- 9 files changed, 21 insertions(+), 21 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp index 95548579..78e00e57 100644 --- a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp @@ -37,11 +37,11 @@ int main() { v2.emplace_back(3.2, std::numeric_limits::infinity()); - double b = Gudhi::bottleneck_distance::compute(v1, v2); + double b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2); std::cout << "Bottleneck distance = " << b << std::endl; - b = Gudhi::bottleneck_distance::compute(v1, v2, 0.1); + b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1); std::cout << "Approx bottleneck distance = " << b << std::endl; diff --git a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp index 4a503e4c..ceedccc5 100644 --- a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp @@ -72,6 +72,6 @@ int main( int argc , char** argv ) { tolerance = atof( argv[3] ); } - double b = Gudhi::bottleneck_distance::compute(diag1, diag2, tolerance); + double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance); std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; } diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index 4a69fb96..0c9e5ea8 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -28,7 +28,7 @@ namespace Gudhi { -namespace bottleneck_distance { +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. @@ -36,7 +36,7 @@ namespace bottleneck_distance { * \ingroup bottleneck_distance */ template -double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=0.) { +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()) @@ -69,7 +69,7 @@ double compute(const Persistence_diagram1 &diag1, const Persistence_diagram2 &di -} // namespace bottleneck_distance +} // namespace persistence_diagram } // namespace Gudhi diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h index 10b7073a..e9f455d7 100644 --- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -27,7 +27,7 @@ namespace Gudhi { -namespace bottleneck_distance { +namespace persistence_diagram { /** \internal \brief Structure representing a graph matching. The graph is a Persistence_diagrams_graph. * @@ -172,7 +172,7 @@ inline void Graph_matching::update(std::vector& path) { } -} // namespace bottleneck_distance +} // namespace persistence_diagram } // namespace Gudhi diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h index 27a59d3f..70342d0c 100644 --- a/src/Bottleneck_distance/include/gudhi/Internal_point.h +++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h @@ -25,7 +25,7 @@ namespace Gudhi { -namespace bottleneck_distance { +namespace persistence_diagram { /** \internal \brief Returns the used index for encoding none of the points */ int null_point_index(); @@ -59,7 +59,7 @@ struct Construct_coord_iterator { { return p.vec+2; } }; -} // namespace bottleneck_distance +} // namespace persistence_diagram } // namespace Gudhi diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h index e6acafc6..792925b7 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -38,7 +38,7 @@ namespace Gudhi { -namespace bottleneck_distance { +namespace persistence_diagram { /** \internal \brief data structure used to find any point (including projections) in V near to a query point from U * (which can be a projection). @@ -162,7 +162,7 @@ inline int Layered_neighbors_finder::vlayers_number() const { return static_cast(neighbors_finder.size()); } -} // namespace bottleneck_distance +} // namespace persistence_diagram } // namespace Gudhi diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h index f16a6c95..2ee55995 100644 --- a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h +++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h @@ -29,7 +29,7 @@ namespace Gudhi { -namespace bottleneck_distance { +namespace persistence_diagram { /** \internal \brief Structure representing an euclidean bipartite graph containing @@ -172,7 +172,7 @@ inline double Persistence_graph::diameter_bound() const { } -} // namespace bottleneck_distance +} // namespace persistence_diagram } // namespace Gudhi diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp index 8bd41104..f51714a2 100644 --- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp +++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp @@ -25,7 +25,7 @@ #include #include -using namespace Gudhi::bottleneck_distance; +using namespace Gudhi::persistence_diagram; double upper_bound = 400.; // any real >0 @@ -34,7 +34,7 @@ int main(){ std::ofstream objetfichier; objetfichier.open("results.csv", std::ios::out); - for(int n = 400; n<=2000; n+=400){ + for(int n = 1000; n<=4000; n+=1000){ std::uniform_real_distribution unif1(0.,upper_bound); std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); std::default_random_engine re; @@ -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.00001); + double b = bottleneck_distance(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 2e67d436..fba1d369 100644 --- a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp +++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp @@ -28,7 +28,7 @@ #include #include -using namespace Gudhi::bottleneck_distance; +using namespace Gudhi::persistence_diagram; int n1 = 81; // a natural number >0 int n2 = 180; // a natural number >0 @@ -161,7 +161,7 @@ BOOST_AUTO_TEST_CASE(global){ if(i%3==0) 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/10000.) <= upper_bound/100. + upper_bound/10000.); - BOOST_CHECK(std::abs(compute(v1, v2) - compute(v1, v2, upper_bound/10000.)) <= upper_bound/10000.); + BOOST_CHECK(bottleneck_distance(v1, v2, 0.) <= upper_bound/100.); + BOOST_CHECK(bottleneck_distance(v1, v2, upper_bound/10000.) <= upper_bound/100. + upper_bound/10000.); + BOOST_CHECK(std::abs(bottleneck_distance(v1, v2, 0.) - bottleneck_distance(v1, v2, upper_bound/10000.)) <= upper_bound/10000.); } -- cgit v1.2.3 From 93c82310194d5fc28979e86287f15dacddbaa5ba Mon Sep 17 00:00:00 2001 From: fgodi Date: Thu, 15 Dec 2016 16:24:35 +0000 Subject: small modifications git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1882 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 11247904405f3da37514e5cc4848c6971edb4254 --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 5 ++--- src/Bottleneck_distance/test/bottleneck_chrono.cpp | 3 ++- 2 files changed, 4 insertions(+), 4 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index 0c9e5ea8..24a31ac1 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -46,12 +46,11 @@ double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_ sd = g.sorted_distances(); long idmin = 0; long idmax = e==0. ? sd.size() - 1 : g.diameter_bound()/e + 1; - // alpha can change the complexity - double alpha = std::pow(idmax, 0.25); + 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) / alpha); + 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) diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp index f51714a2..dcf98460 100644 --- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp +++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp @@ -51,8 +51,9 @@ int main(){ if(i%3==0) v2.emplace_back(std::max(a,b),std::max(a,b)+y); } + double epsilon = 0.0000000000000001; std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); - double b = bottleneck_distance(v1,v2, 0.00001); + double b = bottleneck_distance(v1,v2, epsilon); 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)); -- cgit v1.2.3 From 07db62d323a2dabdc86ad6f1eb8edbc7b1bf26bb Mon Sep 17 00:00:00 2001 From: fgodi Date: Thu, 15 Dec 2016 16:51:19 +0000 Subject: small modifications git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1885 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 111fc4e381223fc2d4b6e6c71b61dd1552c7cc95 --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 2 +- src/Bottleneck_distance/test/bottleneck_chrono.cpp | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index 24a31ac1..bb1d6d38 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -49,7 +49,7 @@ double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_ double alpha = std::pow(g.size(), 1./5.); Graph_matching m(g); Graph_matching biggest_unperfect(g); - while (idmin != idmax) { + while (idmax - idmin > epsilon) { long step = static_cast((idmax - idmin - 1)/alpha); m.set_r(e == 0. ? sd.at(idmin + step) : e*(idmin + step)); while (m.multi_augment()); diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp index dcf98460..7a440aeb 100644 --- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp +++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp @@ -51,6 +51,7 @@ int main(){ if(i%3==0) v2.emplace_back(std::max(a,b),std::max(a,b)+y); } + std::cout << std::numeric_limits::max() / std::numeric_limits::min() << std::endl; double epsilon = 0.0000000000000001; std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); double b = bottleneck_distance(v1,v2, epsilon); -- cgit v1.2.3 From b93ea27ea392f49f85deee23526c9330a716093b Mon Sep 17 00:00:00 2001 From: fgodi Date: Thu, 15 Dec 2016 16:57:53 +0000 Subject: cancel commit of temporary testing code git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1886 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 03f3e6515fd74fd2ae8d1de6b7febb4af7607fb5 --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 2 +- src/Bottleneck_distance/test/bottleneck_chrono.cpp | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index bb1d6d38..24a31ac1 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -49,7 +49,7 @@ double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_ double alpha = std::pow(g.size(), 1./5.); Graph_matching m(g); Graph_matching biggest_unperfect(g); - while (idmax - idmin > epsilon) { + 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()); diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp index 7a440aeb..dcf98460 100644 --- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp +++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp @@ -51,7 +51,6 @@ int main(){ if(i%3==0) v2.emplace_back(std::max(a,b),std::max(a,b)+y); } - std::cout << std::numeric_limits::max() / std::numeric_limits::min() << std::endl; double epsilon = 0.0000000000000001; std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); double b = bottleneck_distance(v1,v2, epsilon); -- cgit v1.2.3 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 +-- src/Bottleneck_distance/test/bottleneck_chrono.cpp | 2 +- 3 files changed, 50 insertions(+), 31 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') 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); } 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::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(); -- cgit v1.2.3 From 66c657a9b422474a1a151fea2ba711b377dce57f Mon Sep 17 00:00:00 2001 From: fgodi Date: Fri, 16 Dec 2016 12:19:22 +0000 Subject: modifs done git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1897 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: f1c2ac9389d5da7556a5c47e90c68572ae97450b --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 16 +++++++++++----- src/Bottleneck_distance/test/bottleneck_chrono.cpp | 2 +- 2 files changed, 12 insertions(+), 6 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index c34ea933..d7e11a05 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -34,6 +34,12 @@ 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.); + if(e < std::numeric_limits::epsilon() * alpha){ + e = std::numeric_limits::epsilon() * alpha; +#ifdef DEBUG_TRACES + std::cout << "Epsilon user given value is less than eps_min. Forced to eps_min by the application" << std::endl; +#endif // DEBUG_TRACES + } Graph_matching m(g); Graph_matching biggest_unperfect(g); while (b_upper_bound - b_lower_bound > 2*e) { @@ -73,18 +79,18 @@ double bottleneck_distance_exact(Persistence_graph& g) { 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. +/** \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, you get an additive e-approximation, which is a lot faster to compute whatever is e. + * Thus, by default, e is a very small positive double, actually the smallest double possible such that the floating-point inaccuracies don't lead to a failure of the algorithm. * * \ingroup bottleneck_distance */ template -double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=0.) { +double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=std::numeric_limits::epsilon()) { 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); + return std::max(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e)); } } // namespace persistence_diagram diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp index b46a19e4..1bff96f4 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 = std::numeric_limits::min(); + double epsilon = std::numeric_limits::epsilon(); 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(); -- cgit v1.2.3 From b6b00e19be498e039ac40368157b3b7e10da4a0b Mon Sep 17 00:00:00 2001 From: fgodi Date: Fri, 16 Dec 2016 15:36:20 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1903 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: ab471ba8c17c17d4e52a5eaccb19ea36c5f278c8 --- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 10 +++------- src/Bottleneck_distance/test/bottleneck_chrono.cpp | 2 +- 2 files changed, 4 insertions(+), 8 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index d7e11a05..42a0d444 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -34,16 +34,12 @@ 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.); - if(e < std::numeric_limits::epsilon() * alpha){ - e = std::numeric_limits::epsilon() * alpha; -#ifdef DEBUG_TRACES - std::cout << "Epsilon user given value is less than eps_min. Forced to eps_min by the application" << std::endl; -#endif // DEBUG_TRACES - } Graph_matching m(g); Graph_matching biggest_unperfect(g); while (b_upper_bound - b_lower_bound > 2*e) { double step = b_lower_bound + (b_upper_bound - b_lower_bound)/alpha; + if(step <= b_lower_bound || step >= b_upper_bound) //Avoid precision problem + break; m.set_r(step); while (m.multi_augment()); //compute a maximum matching (in the graph corresponding to the current r) if (m.perfect()) { @@ -86,7 +82,7 @@ double bottleneck_distance_exact(Persistence_graph& g) { * \ingroup bottleneck_distance */ template -double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=std::numeric_limits::epsilon()) { +double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=std::numeric_limits::min()) { Persistence_graph g(diag1, diag2, e); if(g.bottleneck_alive() == std::numeric_limits::infinity()) return std::numeric_limits::infinity(); diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp index 0a3bb012..a30d42b5 100644 --- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp +++ b/src/Bottleneck_distance/test/bottleneck_chrono.cpp @@ -34,7 +34,7 @@ int main(){ std::ofstream objetfichier; objetfichier.open("results.csv", std::ios::out); - for(int n = 1000; n<=4000; n+=1000){ + for(int n = 1000; n<=10000; n+=1000){ std::uniform_real_distribution unif1(0.,upper_bound); std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); std::default_random_engine re; -- cgit v1.2.3 From 64dd2a2b0eec1374ed23ca079f86b312125d03f7 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Wed, 4 Jan 2017 08:39:18 +0000 Subject: Move bottleneck_chrono in benchmark Add test in cmake for basic example Move CGAL gudhi patches for bottleneck in dedicated directory Fix cpplint syntax issue git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_integration@1920 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 24671facf791de93dc6fd94bb39ca7362bb22959 --- CMakeLists.txt | 1 + src/Bottleneck_distance/benchmark/CMakeLists.txt | 13 + .../benchmark/bottleneck_chrono.cpp | 62 ++ .../doc/Intro_bottleneck_distance.h | 7 +- src/Bottleneck_distance/example/CMakeLists.txt | 2 + .../example/bottleneck_basic_example.cpp | 26 +- .../example/bottleneck_read_file_example.cpp | 71 ++- src/Bottleneck_distance/include/gudhi/Bottleneck.h | 93 +-- .../include/gudhi/CGAL/Kd_tree.h | 582 ------------------- .../include/gudhi/CGAL/Kd_tree_node.h | 586 ------------------- .../CGAL/Orthogonal_incremental_neighbor_search.h | 621 --------------------- .../include/gudhi/Graph_matching.h | 220 ++++---- .../include/gudhi/Internal_point.h | 69 ++- .../include/gudhi/Neighbors_finder.h | 172 +++--- .../include/gudhi/Persistence_graph.h | 191 ++++--- src/Bottleneck_distance/test/CMakeLists.txt | 1 - src/Bottleneck_distance/test/bottleneck_chrono.cpp | 62 -- .../test/bottleneck_unit_test.cpp | 238 ++++---- src/Tangential_complex/benchmark/CMakeLists.txt | 9 - .../Bottleneck_distance_CGAL_patches.txt | 3 + src/common/include/gudhi_patches/CGAL/Kd_tree.h | 582 +++++++++++++++++++ .../include/gudhi_patches/CGAL/Kd_tree_node.h | 586 +++++++++++++++++++ .../CGAL/Orthogonal_incremental_neighbor_search.h | 621 +++++++++++++++++++++ .../Tangential_complex_CGAL_patches.txt | 82 +++ 24 files changed, 2504 insertions(+), 2396 deletions(-) create mode 100644 src/Bottleneck_distance/benchmark/CMakeLists.txt create mode 100644 src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp delete mode 100644 src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h delete mode 100644 src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h delete mode 100644 src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h delete mode 100644 src/Bottleneck_distance/test/bottleneck_chrono.cpp create mode 100644 src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt create mode 100644 src/common/include/gudhi_patches/CGAL/Kd_tree.h create mode 100644 src/common/include/gudhi_patches/CGAL/Kd_tree_node.h create mode 100644 src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h create mode 100644 src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/CMakeLists.txt b/CMakeLists.txt index d1a3c7bf..f949a803 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -158,6 +158,7 @@ else() add_subdirectory(src/Tangential_complex/benchmark) add_subdirectory(src/Bottleneck_distance/example) add_subdirectory(src/Bottleneck_distance/test) + add_subdirectory(src/Bottleneck_distance/benchmark) # data points generator add_subdirectory(data/points/generator) diff --git a/src/Bottleneck_distance/benchmark/CMakeLists.txt b/src/Bottleneck_distance/benchmark/CMakeLists.txt new file mode 100644 index 00000000..f70dd8ff --- /dev/null +++ b/src/Bottleneck_distance/benchmark/CMakeLists.txt @@ -0,0 +1,13 @@ +cmake_minimum_required(VERSION 2.6) +project(Bottleneck_distance_benchmark) + + +# requires CGAL 4.8 +# cmake -DCGAL_DIR=~/workspace/CGAL-4.8 ../../.. +if(CGAL_FOUND) + if (NOT CGAL_VERSION VERSION_LESS 4.8.0) + if (EIGEN3_FOUND) + add_executable ( bottleneck_chrono bottleneck_chrono.cpp ) + endif() + endif () +endif() diff --git a/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp new file mode 100644 index 00000000..456c570b --- /dev/null +++ b/src/Bottleneck_distance/benchmark/bottleneck_chrono.cpp @@ -0,0 +1,62 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: Francois Godi + * + * Copyright (C) 2015 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +#include + +using namespace Gudhi::persistence_diagram; + + +double upper_bound = 400.; // any real > 0 + +int main() { + std::ofstream result_file; + result_file.open("results.csv", std::ios::out); + + for (int n = 1000; n <= 10000; n += 1000) { + std::uniform_real_distribution unif1(0., upper_bound); + std::uniform_real_distribution unif2(upper_bound / 1000., upper_bound / 100.); + std::default_random_engine re; + std::vector< std::pair > v1, v2; + for (int i = 0; i < n; 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); + } + std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); + double b = bottleneck_distance(v1, v2); + 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)); + result_file << n << ";" << duration.count() << ";" << b << std::endl; + } + result_file.close(); +} diff --git a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h index ebe1123b..21187f9c 100644 --- a/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h +++ b/src/Bottleneck_distance/doc/Intro_bottleneck_distance.h @@ -4,7 +4,7 @@ * * Author: François Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -20,8 +20,8 @@ * along with this program. If not, see . */ -#ifndef DOC_BOTTLENECK_DISTANCE_H_ -#define DOC_BOTTLENECK_DISTANCE_H_ +#ifndef DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ +#define DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ // needs namespace for Doxygen to link on classes namespace Gudhi { @@ -48,3 +48,4 @@ namespace bottleneck_distance { } // namespace Gudhi +#endif // DOC_BOTTLENECK_DISTANCE_INTRO_BOTTLENECK_DISTANCE_H_ diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt index cd53ccfc..c66623e9 100644 --- a/src/Bottleneck_distance/example/CMakeLists.txt +++ b/src/Bottleneck_distance/example/CMakeLists.txt @@ -8,6 +8,8 @@ if(CGAL_FOUND) if (EIGEN3_FOUND) add_executable (bottleneck_read_file_example bottleneck_read_file_example.cpp) add_executable (bottleneck_basic_example bottleneck_basic_example.cpp) + + add_test(bottleneck_basic_example ${CMAKE_CURRENT_BINARY_DIR}/bottleneck_basic_example) endif() endif () endif() diff --git a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp index 78e00e57..91a7302f 100644 --- a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp @@ -4,7 +4,7 @@ * * Authors: Francois Godi, small modifications by Pawel Dlotko * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -25,24 +25,24 @@ int main() { - std::vector< std::pair > v1, v2; + std::vector< std::pair > v1, v2; - v1.emplace_back(2.7, 3.7); - v1.emplace_back(9.6, 14.); - v1.emplace_back(34.2, 34.974); - v1.emplace_back(3., std::numeric_limits::infinity()); + v1.emplace_back(2.7, 3.7); + v1.emplace_back(9.6, 14.); + v1.emplace_back(34.2, 34.974); + v1.emplace_back(3., std::numeric_limits::infinity()); - v2.emplace_back(2.8, 4.45); - v2.emplace_back(9.5, 14.1); - v2.emplace_back(3.2, std::numeric_limits::infinity()); + v2.emplace_back(2.8, 4.45); + v2.emplace_back(9.5, 14.1); + v2.emplace_back(3.2, std::numeric_limits::infinity()); - double b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2); + double b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2); - std::cout << "Bottleneck distance = " << b << std::endl; + std::cout << "Bottleneck distance = " << b << std::endl; - b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1); + b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1); - std::cout << "Approx bottleneck distance = " << b << std::endl; + std::cout << "Approx bottleneck distance = " << b << std::endl; } diff --git a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp index ceedccc5..4c74b66e 100644 --- a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp @@ -4,7 +4,7 @@ * * Authors: Francois Godi, small modifications by Pawel Dlotko * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -28,49 +28,42 @@ #include #include -std::vector< std::pair > read_diagram_from_file( const char* filename ) -{ - std::ifstream in; - in.open( filename ); - std::vector< std::pair > result; - if ( !in.is_open() ) - { - std::cerr << "File : " << filename << " do not exist. The program will now terminate \n"; - throw "File do not exist \n"; - } +std::vector< std::pair > read_diagram_from_file(const char* filename) { + std::ifstream in; + in.open(filename); + std::vector< std::pair > result; + if (!in.is_open()) { + std::cerr << "File : " << filename << " do not exist. The program will now terminate \n"; + throw "File do not exist \n"; + } - std::string line; - while (!in.eof()) - { - getline(in,line); - if ( line.length() != 0 ) - { - std::stringstream lineSS; - lineSS << line; - double beginn, endd; - lineSS >> beginn; - lineSS >> endd; - result.push_back( std::make_pair( beginn , endd ) ); - } + std::string line; + while (!in.eof()) { + getline(in, line); + if (line.length() != 0) { + std::stringstream lineSS; + lineSS << line; + double beginn, endd; + lineSS >> beginn; + lineSS >> endd; + result.push_back(std::make_pair(beginn, endd)); } - in.close(); - return result; -} //read_diagram_from_file + } + in.close(); + return result; +} // read_diagram_from_file -int main( int argc , char** argv ) -{ - if ( argc < 3 ) - { - std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file " << - "should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck" << - " distance (set by default to zero). The program will now terminate \n"; +int main(int argc, char** argv) { + if (argc < 3) { + std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file " << + "should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck" << + " distance (set by default to zero). The program will now terminate \n"; } - std::vector< std::pair< double , double > > diag1 = read_diagram_from_file( argv[1] ); - std::vector< std::pair< double , double > > diag2 = read_diagram_from_file( argv[2] ); + std::vector< std::pair< double, double > > diag1 = read_diagram_from_file(argv[1]); + std::vector< std::pair< double, double > > diag2 = read_diagram_from_file(argv[2]); double tolerance = 0.; - if ( argc == 4 ) - { - tolerance = atof( argv[3] ); + if (argc == 4) { + tolerance = atof(argv[3]); } double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance); std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index 42a0d444..2b7e4767 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -31,62 +31,65 @@ namespace Gudhi { namespace persistence_diagram { 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 (b_upper_bound - b_lower_bound > 2*e) { - double step = b_lower_bound + (b_upper_bound - b_lower_bound)/alpha; - if(step <= b_lower_bound || step >= b_upper_bound) //Avoid precision problem - break; - m.set_r(step); - while (m.multi_augment()); //compute a maximum matching (in the graph corresponding to the current r) - if (m.perfect()) { - m = biggest_unperfect; - b_upper_bound = step; - } else { - biggest_unperfect = m; - b_lower_bound = step; - } + 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 (b_upper_bound - b_lower_bound > 2 * e) { + double step = b_lower_bound + (b_upper_bound - b_lower_bound) / alpha; + if (step <= b_lower_bound || step >= b_upper_bound) // Avoid precision problem + break; + m.set_r(step); + while (m.multi_augment()); // compute a maximum matching (in the graph corresponding to the current r) + if (m.perfect()) { + m = biggest_unperfect; + b_upper_bound = step; + } else { + biggest_unperfect = m; + b_lower_bound = step; } - return (b_lower_bound + b_upper_bound)/2.; + } + 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; - } + 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); + } + 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, you get an additive e-approximation, which is a lot faster to compute whatever is e. - * Thus, by default, e is a very small positive double, actually the smallest double possible such that the floating-point inaccuracies don't lead to a failure of the algorithm. + * If the last parameter e is not 0, you get an additive e-approximation, which is a lot faster to compute whatever is + * e. + * Thus, by default, e is a very small positive double, actually the smallest double possible such that the + * floating-point inaccuracies don't lead to a failure of the algorithm. * * \ingroup bottleneck_distance */ template -double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=std::numeric_limits::min()) { - Persistence_graph g(diag1, diag2, e); - if(g.bottleneck_alive() == std::numeric_limits::infinity()) - return std::numeric_limits::infinity(); - return std::max(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e)); +double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, + double e = std::numeric_limits::min()) { + Persistence_graph g(diag1, diag2, e); + if (g.bottleneck_alive() == std::numeric_limits::infinity()) + return std::numeric_limits::infinity(); + return std::max(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e)); } } // namespace persistence_diagram diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h deleted file mode 100644 index f085b0da..00000000 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree.h +++ /dev/null @@ -1,582 +0,0 @@ -// Copyright (c) 2002,2011,2014 Utrecht University (The Netherlands), Max-Planck-Institute Saarbruecken (Germany). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// You can redistribute it and/or modify it under the terms of the GNU -// General Public License as published by the Free Software Foundation, -// either version 3 of the License, or (at your option) any later version. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL$ -// $Id$ -// -// Author(s) : Hans Tangelder (), -// : Waqar Khan - -#ifndef CGAL_KD_TREE_H -#define CGAL_KD_TREE_H - -#include "Kd_tree_node.h" - -#include -#include -#include - -#include -#include -#include - - -#include -#include -#include - -#ifdef CGAL_HAS_THREADS -#include -#endif - -namespace CGAL { - -//template , class UseExtendedNode = Tag_true > -template , class UseExtendedNode = Tag_true > -class Kd_tree { - -public: - typedef SearchTraits Traits; - typedef Splitter_ Splitter; - typedef typename SearchTraits::Point_d Point_d; - typedef typename Splitter::Container Point_container; - - typedef typename SearchTraits::FT FT; - typedef Kd_tree_node Node; - typedef Kd_tree_leaf_node Leaf_node; - typedef Kd_tree_internal_node Internal_node; - typedef Kd_tree Tree; - typedef Kd_tree Self; - - typedef Node* Node_handle; - typedef const Node* Node_const_handle; - typedef Leaf_node* Leaf_node_handle; - typedef const Leaf_node* Leaf_node_const_handle; - typedef Internal_node* Internal_node_handle; - typedef const Internal_node* Internal_node_const_handle; - typedef typename std::vector::const_iterator Point_d_iterator; - typedef typename std::vector::const_iterator Point_d_const_iterator; - typedef typename Splitter::Separator Separator; - typedef typename std::vector::const_iterator iterator; - typedef typename std::vector::const_iterator const_iterator; - - typedef typename std::vector::size_type size_type; - - typedef typename internal::Get_dimension_tag::Dimension D; - -private: - SearchTraits traits_; - Splitter split; - - - // wokaround for https://svn.boost.org/trac/boost/ticket/9332 -#if (_MSC_VER == 1800) && (BOOST_VERSION == 105500) - std::deque internal_nodes; - std::deque leaf_nodes; -#else - boost::container::deque internal_nodes; - boost::container::deque leaf_nodes; -#endif - - Node_handle tree_root; - - Kd_tree_rectangle* bbox; - std::vector pts; - - // Instead of storing the points in arrays in the Kd_tree_node - // we put all the data in a vector in the Kd_tree. - // and we only store an iterator range in the Kd_tree_node. - // - std::vector data; - - - #ifdef CGAL_HAS_THREADS - mutable CGAL_MUTEX building_mutex;//mutex used to protect const calls inducing build() - #endif - bool built_; - bool removed_; - - // protected copy constructor - Kd_tree(const Tree& tree) - : traits_(tree.traits_),built_(tree.built_) - {}; - - - // Instead of the recursive construction of the tree in the class Kd_tree_node - // we do this in the tree class. The advantage is that we then can optimize - // the allocation of the nodes. - - // The leaf node - Node_handle - create_leaf_node(Point_container& c) - { - Leaf_node node(true , static_cast(c.size())); - std::ptrdiff_t tmp = c.begin() - data.begin(); - node.data = pts.begin() + tmp; - - leaf_nodes.push_back(node); - Leaf_node_handle nh = &leaf_nodes.back(); - - - return nh; - } - - - // The internal node - - Node_handle - create_internal_node(Point_container& c, const Tag_true&) - { - return create_internal_node_use_extension(c); - } - - Node_handle - create_internal_node(Point_container& c, const Tag_false&) - { - return create_internal_node(c); - } - - - - // TODO: Similiar to the leaf_init function above, a part of the code should be - // moved to a the class Kd_tree_node. - // It is not proper yet, but the goal was to see if there is - // a potential performance gain through the Compact_container - Node_handle - create_internal_node_use_extension(Point_container& c) - { - Internal_node node(false); - internal_nodes.push_back(node); - Internal_node_handle nh = &internal_nodes.back(); - - Separator sep; - Point_container c_low(c.dimension(),traits_); - split(sep, c, c_low); - nh->set_separator(sep); - - int cd = nh->cutting_dimension(); - if(!c_low.empty()){ - nh->lower_low_val = c_low.tight_bounding_box().min_coord(cd); - nh->lower_high_val = c_low.tight_bounding_box().max_coord(cd); - } - else{ - nh->lower_low_val = nh->cutting_value(); - nh->lower_high_val = nh->cutting_value(); - } - if(!c.empty()){ - nh->upper_low_val = c.tight_bounding_box().min_coord(cd); - nh->upper_high_val = c.tight_bounding_box().max_coord(cd); - } - else{ - nh->upper_low_val = nh->cutting_value(); - nh->upper_high_val = nh->cutting_value(); - } - - CGAL_assertion(nh->cutting_value() >= nh->lower_low_val); - CGAL_assertion(nh->cutting_value() <= nh->upper_high_val); - - if (c_low.size() > split.bucket_size()){ - nh->lower_ch = create_internal_node_use_extension(c_low); - }else{ - nh->lower_ch = create_leaf_node(c_low); - } - if (c.size() > split.bucket_size()){ - nh->upper_ch = create_internal_node_use_extension(c); - }else{ - nh->upper_ch = create_leaf_node(c); - } - - - - - return nh; - } - - - // Note also that I duplicated the code to get rid if the if's for - // the boolean use_extension which was constant over the construction - Node_handle - create_internal_node(Point_container& c) - { - Internal_node node(false); - internal_nodes.push_back(node); - Internal_node_handle nh = &internal_nodes.back(); - Separator sep; - - Point_container c_low(c.dimension(),traits_); - split(sep, c, c_low); - nh->set_separator(sep); - - if (c_low.size() > split.bucket_size()){ - nh->lower_ch = create_internal_node(c_low); - }else{ - nh->lower_ch = create_leaf_node(c_low); - } - if (c.size() > split.bucket_size()){ - nh->upper_ch = create_internal_node(c); - }else{ - nh->upper_ch = create_leaf_node(c); - } - - - - return nh; - } - - - -public: - - Kd_tree(Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) - : traits_(traits),split(s), built_(false), removed_(false) - {} - - template - Kd_tree(InputIterator first, InputIterator beyond, - Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) - : traits_(traits),split(s), built_(false), removed_(false) - { - pts.insert(pts.end(), first, beyond); - } - - bool empty() const { - return pts.empty(); - } - - void - build() - { - // This function is not ready to be called when a tree already exists, one - // must call invalidate_built() first. - CGAL_assertion(!is_built()); - CGAL_assertion(!removed_); - const Point_d& p = *pts.begin(); - typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits_.construct_cartesian_const_iterator_d_object(); - int dim = static_cast(std::distance(ccci(p), ccci(p,0))); - - data.reserve(pts.size()); - for(unsigned int i = 0; i < pts.size(); i++){ - data.push_back(&pts[i]); - } - Point_container c(dim, data.begin(), data.end(),traits_); - bbox = new Kd_tree_rectangle(c.bounding_box()); - if (c.size() <= split.bucket_size()){ - tree_root = create_leaf_node(c); - }else { - tree_root = create_internal_node(c, UseExtendedNode()); - } - - //Reorder vector for spatial locality - std::vector ptstmp; - ptstmp.resize(pts.size()); - for (std::size_t i = 0; i < pts.size(); ++i){ - ptstmp[i] = *data[i]; - } - for(std::size_t i = 0; i < leaf_nodes.size(); ++i){ - std::ptrdiff_t tmp = leaf_nodes[i].begin() - pts.begin(); - leaf_nodes[i].data = ptstmp.begin() + tmp; - } - pts.swap(ptstmp); - - data.clear(); - - built_ = true; - } - -private: - //any call to this function is for the moment not threadsafe - void const_build() const { - #ifdef CGAL_HAS_THREADS - //this ensure that build() will be called once - CGAL_SCOPED_LOCK(building_mutex); - if(!is_built()) - #endif - const_cast(this)->build(); //THIS IS NOT THREADSAFE - } -public: - - bool is_built() const - { - return built_; - } - - void invalidate_built() - { - if(removed_){ - // Walk the tree to collect the remaining points. - // Writing directly to pts would likely work, but better be safe. - std::vector ptstmp; - //ptstmp.resize(root()->num_items()); - root()->tree_items(std::back_inserter(ptstmp)); - pts.swap(ptstmp); - removed_=false; - CGAL_assertion(is_built()); // the rest of the cleanup must happen - } - if(is_built()){ - internal_nodes.clear(); - leaf_nodes.clear(); - data.clear(); - delete bbox; - built_ = false; - } - } - - void clear() - { - invalidate_built(); - pts.clear(); - removed_ = false; - } - - void - insert(const Point_d& p) - { - invalidate_built(); - pts.push_back(p); - } - - template - void - insert(InputIterator first, InputIterator beyond) - { - invalidate_built(); - pts.insert(pts.end(),first, beyond); - } - -private: - struct Equal_by_coordinates { - SearchTraits const* traits; - Point_d const* pp; - bool operator()(Point_d const&q) const { - typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits->construct_cartesian_const_iterator_d_object(); - return std::equal(ccci(*pp), ccci(*pp,0), ccci(q)); - } - }; - Equal_by_coordinates equal_by_coordinates(Point_d const&p){ - Equal_by_coordinates ret = { &traits(), &p }; - return ret; - } - -public: - void - remove(const Point_d& p) - { - remove(p, equal_by_coordinates(p)); - } - - template - void - remove(const Point_d& p, Equal const& equal_to_p) - { -#if 0 - // This code could have quadratic runtime. - if (!is_built()) { - std::vector::iterator pi = std::find(pts.begin(), pts.end(), p); - // Precondition: the point must be there. - CGAL_assertion (pi != pts.end()); - pts.erase(pi); - return; - } -#endif - bool success = remove_(p, 0, false, 0, false, root(), equal_to_p); - CGAL_assertion(success); - - // Do not set the flag is the tree has been cleared. - if(is_built()) - removed_ |= success; - } -private: - template - bool remove_(const Point_d& p, - Internal_node_handle grandparent, bool parent_islower, - Internal_node_handle parent, bool islower, - Node_handle node, Equal const& equal_to_p) { - // Recurse to locate the point - if (!node->is_leaf()) { - Internal_node_handle newparent = static_cast(node); - // FIXME: This should be if(xcutting_dimension()] <= newparent->cutting_value()) { - if (remove_(p, parent, islower, newparent, true, newparent->lower(), equal_to_p)) - return true; - } - //if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] >= newparent->cutting_value()) - return remove_(p, parent, islower, newparent, false, newparent->upper(), equal_to_p); - - CGAL_assertion(false); // Point was not found - } - - // Actual removal - Leaf_node_handle lnode = static_cast(node); - if (lnode->size() > 1) { - iterator pi = std::find_if(lnode->begin(), lnode->end(), equal_to_p); - // FIXME: we should ensure this never happens - if (pi == lnode->end()) return false; - iterator lasti = lnode->end() - 1; - if (pi != lasti) { - // Hack to get a non-const iterator - std::iter_swap(pts.begin()+(pi-pts.begin()), pts.begin()+(lasti-pts.begin())); - } - lnode->drop_last_point(); - } else if (!equal_to_p(*lnode->begin())) { - // FIXME: we should ensure this never happens - return false; - } else if (grandparent) { - Node_handle brother = islower ? parent->upper() : parent->lower(); - if (parent_islower) - grandparent->set_lower(brother); - else - grandparent->set_upper(brother); - } else if (parent) { - tree_root = islower ? parent->upper() : parent->lower(); - } else { - clear(); - } - return true; - } - -public: - //For efficiency; reserve the size of the points vectors in advance (if the number of points is already known). - void reserve(size_t size) - { - pts.reserve(size); - } - - //Get the capacity of the underlying points vector. - size_t capacity() - { - return pts.capacity(); - } - - - template - OutputIterator - search(OutputIterator it, const FuzzyQueryItem& q) const - { - if(! pts.empty()){ - - if(! is_built()){ - const_build(); - } - Kd_tree_rectangle b(*bbox); - return tree_root->search(it,q,b); - } - return it; - } - - - template - boost::optional - search_any_point(const FuzzyQueryItem& q) const - { - if(! pts.empty()){ - - if(! is_built()){ - const_build(); - } - Kd_tree_rectangle b(*bbox); - return tree_root->search_any_point(q,b); - } - return boost::none; - } - - - ~Kd_tree() { - if(is_built()){ - delete bbox; - } - } - - - const SearchTraits& - traits() const - { - return traits_; - } - - Node_const_handle - root() const - { - if(! is_built()){ - const_build(); - } - return tree_root; - } - - Node_handle - root() - { - if(! is_built()){ - build(); - } - return tree_root; - } - - void - print() const - { - if(! is_built()){ - const_build(); - } - root()->print(); - } - - const Kd_tree_rectangle& - bounding_box() const - { - if(! is_built()){ - const_build(); - } - return *bbox; - } - - const_iterator - begin() const - { - return pts.begin(); - } - - const_iterator - end() const - { - return pts.end(); - } - - size_type - size() const - { - return pts.size(); - } - - // Print statistics of the tree. - std::ostream& - statistics(std::ostream& s) const - { - if(! is_built()){ - const_build(); - } - s << "Tree statistics:" << std::endl; - s << "Number of items stored: " - << root()->num_items() << std::endl; - s << "Number of nodes: " - << root()->num_nodes() << std::endl; - s << " Tree depth: " << root()->depth() << std::endl; - return s; - } - - -}; - -} // namespace CGAL - -#endif // CGAL_KD_TREE_H diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h b/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h deleted file mode 100644 index 909ee260..00000000 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Kd_tree_node.h +++ /dev/null @@ -1,586 +0,0 @@ -// Copyright (c) 2002,2011 Utrecht University (The Netherlands). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// You can redistribute it and/or modify it under the terms of the GNU -// General Public License as published by the Free Software Foundation, -// either version 3 of the License, or (at your option) any later version. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL$ -// $Id$ -// -// -// Authors : Hans Tangelder () - -#ifndef CGAL_KD_TREE_NODE_H -#define CGAL_KD_TREE_NODE_H - -#include "CGAL/Splitters.h" - -#include -#include - -namespace CGAL { - - template - class Kd_tree; - - template < class TreeTraits, class Splitter, class UseExtendedNode > - class Kd_tree_node { - - friend class Kd_tree; - - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; - typedef typename Kd_tree::Internal_node_handle Internal_node_handle; - typedef typename Kd_tree::Internal_node_const_handle Internal_node_const_handle; - typedef typename Kd_tree::Leaf_node_handle Leaf_node_handle; - typedef typename Kd_tree::Leaf_node_const_handle Leaf_node_const_handle; - typedef typename TreeTraits::Point_d Point_d; - - typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::Separator Separator; - typedef typename Kd_tree::Point_d_iterator Point_d_iterator; - typedef typename Kd_tree::iterator iterator; - typedef typename Kd_tree::D D; - - bool leaf; - - public : - Kd_tree_node(bool leaf_) - :leaf(leaf_){} - - bool is_leaf() const{ - return leaf; - } - - std::size_t - num_items() const - { - if (is_leaf()){ - Leaf_node_const_handle node = - static_cast(this); - return node->size(); - } - else { - Internal_node_const_handle node = - static_cast(this); - return node->lower()->num_items() + node->upper()->num_items(); - } - } - - std::size_t - num_nodes() const - { - if (is_leaf()) return 1; - else { - Internal_node_const_handle node = - static_cast(this); - return node->lower()->num_nodes() + node->upper()->num_nodes(); - } - } - - int - depth(const int current_max_depth) const - { - if (is_leaf()){ - return current_max_depth; - } - else { - Internal_node_const_handle node = - static_cast(this); - return - (std::max)( node->lower()->depth(current_max_depth + 1), - node->upper()->depth(current_max_depth + 1)); - } - } - - int - depth() const - { - return depth(1); - } - - template - OutputIterator - tree_items(OutputIterator it) const { - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(this); - if (node->size()>0) - for (iterator i=node->begin(); i != node->end(); i++) - {*it=*i; ++it;} - } - else { - Internal_node_const_handle node = - static_cast(this); - it=node->lower()->tree_items(it); - it=node->upper()->tree_items(it); - } - return it; - } - - - boost::optional - any_tree_item() const { - boost::optional result = boost::none; - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(this); - if (node->size()>0){ - return boost::make_optional(*(node->begin())); - } - } - else { - Internal_node_const_handle node = - static_cast(this); - result = node->lower()->any_tree_item(); - if(! result){ - result = node->upper()->any_tree_item(); - } - } - return result; - } - - - void - indent(int d) const - { - for(int i = 0; i < d; i++){ - std::cout << " "; - } - } - - - void - print(int d = 0) const - { - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(this); - indent(d); - std::cout << "leaf" << std::endl; - if (node->size()>0) - for (iterator i=node->begin(); i != node->end(); i++) - {indent(d);std::cout << *i << std::endl;} - } - else { - Internal_node_const_handle node = - static_cast(this); - indent(d); - std::cout << "lower tree" << std::endl; - node->lower()->print(d+1); - indent(d); - std::cout << "separator: dim = " << node->cutting_dimension() << " val = " << node->cutting_value() << std::endl; - indent(d); - std::cout << "upper tree" << std::endl; - node->upper()->print(d+1); - } - } - - - template - OutputIterator - search(OutputIterator it, const FuzzyQueryItem& q, - Kd_tree_rectangle& b) const - { - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(this); - if (node->size()>0) - for (iterator i=node->begin(); i != node->end(); i++) - if (q.contains(*i)) - {*it++=*i;} - } - else { - Internal_node_const_handle node = - static_cast(this); - // after splitting b denotes the lower part of b - Kd_tree_rectangle b_upper(b); - b.split(b_upper, node->cutting_dimension(), - node->cutting_value()); - - if (q.outer_range_contains(b)) - it=node->lower()->tree_items(it); - else - if (q.inner_range_intersects(b)) - it=node->lower()->search(it,q,b); - if (q.outer_range_contains(b_upper)) - it=node->upper()->tree_items(it); - else - if (q.inner_range_intersects(b_upper)) - it=node->upper()->search(it,q,b_upper); - }; - return it; - } - - - template - boost::optional - search_any_point(const FuzzyQueryItem& q, - Kd_tree_rectangle& b) const - { - boost::optional result = boost::none; - if (is_leaf()) { - Leaf_node_const_handle node = - static_cast(this); - if (node->size()>0) - for (iterator i=node->begin(); i != node->end(); i++) - if (q.contains(*i)) - { result = *i; break; } - } - else { - Internal_node_const_handle node = - static_cast(this); - // after splitting b denotes the lower part of b - Kd_tree_rectangle b_upper(b); - b.split(b_upper, node->cutting_dimension(), - node->cutting_value()); - - if (q.outer_range_contains(b)){ - result = node->lower()->any_tree_item(); - }else{ - if (q.inner_range_intersects(b)){ - result = node->lower()->search_any_point(q,b); - } - } - if(result){ - return result; - } - if (q.outer_range_contains(b_upper)){ - result = node->upper()->any_tree_item(); - }else{ - if (q.inner_range_intersects(b_upper)) - result = node->upper()->search_any_point(q,b_upper); - } - } - return result; - } - - }; - - - template < class TreeTraits, class Splitter, class UseExtendedNode > - class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ - - friend class Kd_tree; - - typedef typename Kd_tree::iterator iterator; - typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; - typedef typename TreeTraits::Point_d Point_d; - - private: - - // private variables for leaf nodes - boost::int32_t n; // denotes number of items in a leaf node - iterator data; // iterator to data in leaf node - - - public: - - // default constructor - Kd_tree_leaf_node() - {} - - Kd_tree_leaf_node(bool leaf_ ) - : Base(leaf_) - {} - - Kd_tree_leaf_node(bool leaf_,unsigned int n_ ) - : Base(leaf_), n(n_) - {} - - // members for all nodes - - // members for leaf nodes only - inline - unsigned int - size() const - { - return n; - } - - inline - iterator - begin() const - { - return data; - } - - inline - iterator - end() const - { - return data + n; - } - - inline - void - drop_last_point() - { - --n; - } - - }; //leaf node - - - - template < class TreeTraits, class Splitter, class UseExtendedNode> - class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ - - friend class Kd_tree; - - typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; - - typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::Separator Separator; - - private: - - // private variables for internal nodes - boost::int32_t cut_dim; - FT cut_val; - Node_handle lower_ch, upper_ch; - - - // private variables for extended internal nodes - FT upper_low_val; - FT upper_high_val; - FT lower_low_val; - FT lower_high_val; - - - public: - - // default constructor - Kd_tree_internal_node() - {} - - Kd_tree_internal_node(bool leaf_) - : Base(leaf_) - {} - - - // members for internal node and extended internal node - - inline - Node_const_handle - lower() const - { - return lower_ch; - } - - inline - Node_const_handle - upper() const - { - return upper_ch; - } - - inline - Node_handle - lower() - { - return lower_ch; - } - - inline - Node_handle - upper() - { - return upper_ch; - } - - inline - void - set_lower(Node_handle nh) - { - lower_ch = nh; - } - - inline - void - set_upper(Node_handle nh) - { - upper_ch = nh; - } - - // inline Separator& separator() {return sep; } - // use instead - inline - void set_separator(Separator& sep){ - cut_dim = sep.cutting_dimension(); - cut_val = sep.cutting_value(); - } - - inline - FT - cutting_value() const - { - return cut_val; - } - - inline - int - cutting_dimension() const - { - return cut_dim; - } - - // members for extended internal node only - inline - FT - upper_low_value() const - { - return upper_low_val; - } - - inline - FT - upper_high_value() const - { - return upper_high_val; - } - - inline - FT - lower_low_value() const - { - return lower_low_val; - } - - inline - FT - lower_high_value() const - { - return lower_high_val; - } - - /*Separator& - separator() - { - return Separator(cutting_dimension,cutting_value); - }*/ - - - };//internal node - - template < class TreeTraits, class Splitter> - class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, Tag_false >{ - - friend class Kd_tree; - - typedef Kd_tree_node< TreeTraits, Splitter, Tag_false> Base; - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; - - typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::Separator Separator; - - private: - - // private variables for internal nodes - boost::uint8_t cut_dim; - FT cut_val; - - Node_handle lower_ch, upper_ch; - - public: - - // default constructor - Kd_tree_internal_node() - {} - - Kd_tree_internal_node(bool leaf_) - : Base(leaf_) - {} - - - // members for internal node and extended internal node - - inline - Node_const_handle - lower() const - { - return lower_ch; - } - - inline - Node_const_handle - upper() const - { - return upper_ch; - } - - inline - Node_handle - lower() - { - return lower_ch; - } - - inline - Node_handle - upper() - { - return upper_ch; - } - - inline - void - set_lower(Node_handle nh) - { - lower_ch = nh; - } - - inline - void - set_upper(Node_handle nh) - { - upper_ch = nh; - } - - // inline Separator& separator() {return sep; } - // use instead - - inline - void set_separator(Separator& sep){ - cut_dim = sep.cutting_dimension(); - cut_val = sep.cutting_value(); - } - - inline - FT - cutting_value() const - { - return cut_val; - } - - inline - int - cutting_dimension() const - { - return cut_dim; - } - - /* Separator& - separator() - { - return Separator(cutting_dimension,cutting_value); - }*/ - - - };//internal node - - - -} // namespace CGAL -#endif // CGAL_KDTREE_NODE_H diff --git a/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h b/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h deleted file mode 100644 index dbe707ed..00000000 --- a/src/Bottleneck_distance/include/gudhi/CGAL/Orthogonal_incremental_neighbor_search.h +++ /dev/null @@ -1,621 +0,0 @@ -// Copyright (c) 2002,2011 Utrecht University (The Netherlands). -// All rights reserved. -// -// This file is part of CGAL (www.cgal.org). -// You can redistribute it and/or modify it under the terms of the GNU -// General Public License as published by the Free Software Foundation, -// either version 3 of the License, or (at your option) any later version. -// -// Licensees holding a valid commercial license may use this file in -// accordance with the commercial license agreement provided with the software. -// -// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE -// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. -// -// $URL$ -// $Id$ -// -// -// Author(s) : Hans Tangelder () - -#ifndef CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH -#define CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH - -#include "CGAL/Kd_tree.h" - -#include -#include -#include -#include -#include -#include - -namespace CGAL { - - template ::type, - class Splitter_ = Sliding_midpoint, - class Tree_= Kd_tree > - class Orthogonal_incremental_neighbor_search { - - public: - typedef Splitter_ Splitter; - typedef Tree_ Tree; - typedef Distance_ Distance; - typedef typename SearchTraits::Point_d Point_d; - typedef typename Distance::Query_item Query_item; - typedef typename SearchTraits::FT FT; - typedef typename Tree::Point_d_iterator Point_d_iterator; - typedef typename Tree::Node_const_handle Node_const_handle; - - typedef std::pair Point_with_transformed_distance; - typedef CGAL::cpp11::tuple > Node_with_distance; - typedef std::vector Node_with_distance_vector; - typedef std::vector Point_with_transformed_distance_vector; - - template - struct Object_wrapper - { - T object; - Object_wrapper(const T& t):object(t){} - const T& operator* () const { return object; } - const T* operator-> () const { return &object; } - }; - - class Iterator_implementation { - SearchTraits traits; - public: - - int number_of_neighbours_computed; - int number_of_internal_nodes_visited; - int number_of_leaf_nodes_visited; - int number_of_items_visited; - - private: - - typedef std::vector Distance_vector; - - Distance_vector dists; - - Distance Orthogonal_distance_instance; - - FT multiplication_factor; - - Query_item query_point; - - FT distance_to_root; - - bool search_nearest_neighbour; - - FT rd; - - - class Priority_higher { - public: - - bool search_nearest; - - Priority_higher(bool search_the_nearest_neighbour) - : search_nearest(search_the_nearest_neighbour) - {} - - //highest priority is smallest distance - bool - operator() (Node_with_distance* n1, Node_with_distance* n2) const - { - return (search_nearest) ? (CGAL::cpp11::get<1>(*n1) > CGAL::cpp11::get<1>(*n2)) : (CGAL::cpp11::get<1>(*n2) > CGAL::cpp11::get<1>(*n1)); - } - }; - - class Distance_smaller { - - public: - - bool search_nearest; - - Distance_smaller(bool search_the_nearest_neighbour) - : search_nearest(search_the_nearest_neighbour) - {} - - //highest priority is smallest distance - bool operator() (Point_with_transformed_distance* p1, Point_with_transformed_distance* p2) const - { - return (search_nearest) ? (p1->second > p2->second) : (p2->second > p1->second); - } - }; - - - std::priority_queue PriorityQueue; - - public: - std::priority_queue Item_PriorityQueue; - - - public: - - int reference_count; - - - - // constructor - Iterator_implementation(const Tree& tree,const Query_item& q, const Distance& tr, - FT Eps=FT(0.0), bool search_nearest=true) - : traits(tree.traits()),number_of_neighbours_computed(0), number_of_internal_nodes_visited(0), - number_of_leaf_nodes_visited(0), number_of_items_visited(0), - Orthogonal_distance_instance(tr), multiplication_factor(Orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)), - query_point(q), search_nearest_neighbour(search_nearest), - PriorityQueue(Priority_higher(search_nearest)), Item_PriorityQueue(Distance_smaller(search_nearest)), - reference_count(1) - - - { - if (tree.empty()) return; - - typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits.construct_cartesian_const_iterator_d_object(); - int dim = static_cast(std::distance(ccci(q), ccci(q,0))); - - dists.resize(dim); - for(int i=0 ; i(*The_Root); - Compute_the_next_nearest_neighbour(); - } - else{ - distance_to_root= - Orthogonal_distance_instance.max_distance_to_rectangle(q, - tree.bounding_box(), dists); - Node_with_distance *The_Root = new Node_with_distance(tree.root(), - distance_to_root, dists); - PriorityQueue.push(The_Root); - - // rd is the distance of the top of the priority queue to q - rd=CGAL::cpp11::get<1>(*The_Root); - Compute_the_next_furthest_neighbour(); - } - - - } - - // * operator - const Point_with_transformed_distance& - operator* () const - { - return *(Item_PriorityQueue.top()); - } - - // prefix operator - Iterator_implementation& - operator++() - { - Delete_the_current_item_top(); - if(search_nearest_neighbour) - Compute_the_next_nearest_neighbour(); - else - Compute_the_next_furthest_neighbour(); - return *this; - } - - // postfix operator - Object_wrapper - operator++(int) - { - Object_wrapper result( *(Item_PriorityQueue.top()) ); - ++*this; - return result; - } - - // Print statistics of the general priority search process. - std::ostream& - statistics (std::ostream& s) const { - s << "Orthogonal priority search statistics:" - << std::endl; - s << "Number of internal nodes visited:" - << number_of_internal_nodes_visited << std::endl; - s << "Number of leaf nodes visited:" - << number_of_leaf_nodes_visited << std::endl; - s << "Number of items visited:" - << number_of_items_visited << std::endl; - s << "Number of neighbours computed:" - << number_of_neighbours_computed << std::endl; - return s; - } - - - //destructor - ~Iterator_implementation() - { - while (!PriorityQueue.empty()) { - Node_with_distance* The_top=PriorityQueue.top(); - PriorityQueue.pop(); - delete The_top; - } - while (!Item_PriorityQueue.empty()) { - Point_with_transformed_distance* The_top=Item_PriorityQueue.top(); - Item_PriorityQueue.pop(); - delete The_top; - } - } - - private: - - void - Delete_the_current_item_top() - { - Point_with_transformed_distance* The_item_top=Item_PriorityQueue.top(); - Item_PriorityQueue.pop(); - delete The_item_top; - } - - void - Compute_the_next_nearest_neighbour() - { - // compute the next item - bool next_neighbour_found=false; - if (!(Item_PriorityQueue.empty())) { - next_neighbour_found= - (multiplication_factor*rd > Item_PriorityQueue.top()->second); - } - typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); - typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); - // otherwise browse the tree further - while ((!next_neighbour_found) && (!PriorityQueue.empty())) { - Node_with_distance* The_node_top=PriorityQueue.top(); - Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); - dists = CGAL::cpp11::get<2>(*The_node_top); - PriorityQueue.pop(); - delete The_node_top; - FT copy_rd=rd; - while (!(N->is_leaf())) { // compute new distance - typename Tree::Internal_node_const_handle node = - static_cast(N); - number_of_internal_nodes_visited++; - int new_cut_dim=node->cutting_dimension(); - FT new_rd,dst = dists[new_cut_dim]; - FT val = *(query_point_it + new_cut_dim); - FT diff1 = val - node->upper_low_value(); - FT diff2 = val - node->lower_high_value(); - if (diff1 + diff2 < FT(0.0)) { - new_rd= - Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); - - CGAL_assertion(new_rd >= copy_rd); - dists[new_cut_dim] = diff1; - Node_with_distance *Upper_Child = - new Node_with_distance(node->upper(), new_rd, dists); - PriorityQueue.push(Upper_Child); - dists[new_cut_dim] = dst; - N=node->lower(); - - } - else { // compute new distance - new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); - CGAL_assertion(new_rd >= copy_rd); - dists[new_cut_dim] = diff2; - Node_with_distance *Lower_Child = - new Node_with_distance(node->lower(), new_rd, dists); - PriorityQueue.push(Lower_Child); - dists[new_cut_dim] = dst; - N=node->upper(); - } - } - // n is a leaf - typename Tree::Leaf_node_const_handle node = - static_cast(N); - number_of_leaf_nodes_visited++; - if (node->size() > 0) { - for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { - number_of_items_visited++; - FT distance_to_query_point= - Orthogonal_distance_instance.transformed_distance(query_point,*it); - Point_with_transformed_distance *NN_Candidate= - new Point_with_transformed_distance(*it,distance_to_query_point); - Item_PriorityQueue.push(NN_Candidate); - } - // old top of PriorityQueue has been processed, - // hence update rd - - if (!(PriorityQueue.empty())) { - rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); - next_neighbour_found = - (multiplication_factor*rd > - Item_PriorityQueue.top()->second); - } - else // priority queue empty => last neighbour found - { - next_neighbour_found=true; - } - - number_of_neighbours_computed++; - } - } // next_neighbour_found or priority queue is empty - // in the latter case also the item priority quee is empty - } - - - void - Compute_the_next_furthest_neighbour() - { - // compute the next item - bool next_neighbour_found=false; - if (!(Item_PriorityQueue.empty())) { - next_neighbour_found= - (rd < multiplication_factor*Item_PriorityQueue.top()->second); - } - typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); - typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); - // otherwise browse the tree further - while ((!next_neighbour_found) && (!PriorityQueue.empty())) { - Node_with_distance* The_node_top=PriorityQueue.top(); - Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); - dists = CGAL::cpp11::get<2>(*The_node_top); - PriorityQueue.pop(); - delete The_node_top; - FT copy_rd=rd; - while (!(N->is_leaf())) { // compute new distance - typename Tree::Internal_node_const_handle node = - static_cast(N); - number_of_internal_nodes_visited++; - int new_cut_dim=node->cutting_dimension(); - FT new_rd,dst = dists[new_cut_dim]; - FT val = *(query_point_it + new_cut_dim); - FT diff1 = val - node->upper_low_value(); - FT diff2 = val - node->lower_high_value(); - if (diff1 + diff2 < FT(0.0)) { - diff1 = val - node->upper_high_value(); - new_rd= - Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); - Node_with_distance *Lower_Child = - new Node_with_distance(node->lower(), copy_rd, dists); - PriorityQueue.push(Lower_Child); - N=node->upper(); - dists[new_cut_dim] = diff1; - copy_rd=new_rd; - - } - else { // compute new distance - diff2 = val - node->lower_low_value(); - new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); - Node_with_distance *Upper_Child = - new Node_with_distance(node->upper(), copy_rd, dists); - PriorityQueue.push(Upper_Child); - N=node->lower(); - dists[new_cut_dim] = diff2; - copy_rd=new_rd; - } - } - // n is a leaf - typename Tree::Leaf_node_const_handle node = - static_cast(N); - number_of_leaf_nodes_visited++; - if (node->size() > 0) { - for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { - number_of_items_visited++; - FT distance_to_query_point= - Orthogonal_distance_instance.transformed_distance(query_point,*it); - Point_with_transformed_distance *NN_Candidate= - new Point_with_transformed_distance(*it,distance_to_query_point); - Item_PriorityQueue.push(NN_Candidate); - } - // old top of PriorityQueue has been processed, - // hence update rd - - if (!(PriorityQueue.empty())) { - rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); - next_neighbour_found = - (multiplication_factor*rd < - Item_PriorityQueue.top()->second); - } - else // priority queue empty => last neighbour found - { - next_neighbour_found=true; - } - - number_of_neighbours_computed++; - } - } // next_neighbour_found or priority queue is empty - // in the latter case also the item priority quee is empty - } - }; // class Iterator_implementaion - - - - - - - - - - public: - class iterator; - typedef iterator const_iterator; - - // constructor - Orthogonal_incremental_neighbor_search(const Tree& tree, - const Query_item& q, FT Eps = FT(0.0), - bool search_nearest=true, const Distance& tr=Distance()) - : m_tree(tree),m_query(q),m_dist(tr),m_Eps(Eps),m_search_nearest(search_nearest) - {} - - iterator - begin() const - { - return iterator(m_tree,m_query,m_dist,m_Eps,m_search_nearest); - } - - iterator - end() const - { - return iterator(); - } - - std::ostream& - statistics(std::ostream& s) - { - begin()->statistics(s); - return s; - } - - - - - class iterator { - - public: - - typedef std::input_iterator_tag iterator_category; - typedef Point_with_transformed_distance value_type; - typedef Point_with_transformed_distance* pointer; - typedef const Point_with_transformed_distance& reference; - typedef std::size_t size_type; - typedef std::ptrdiff_t difference_type; - typedef int distance_type; - - //class Iterator_implementation; - Iterator_implementation *Ptr_implementation; - - - public: - - // default constructor - iterator() - : Ptr_implementation(0) - {} - - int - the_number_of_items_visited() - { - return Ptr_implementation->number_of_items_visited; - } - - // constructor - iterator(const Tree& tree,const Query_item& q, const Distance& tr=Distance(), FT eps=FT(0.0), - bool search_nearest=true) - : Ptr_implementation(new Iterator_implementation(tree, q, tr, eps, search_nearest)) - {} - - // copy constructor - iterator(const iterator& Iter) - { - Ptr_implementation = Iter.Ptr_implementation; - if (Ptr_implementation != 0) Ptr_implementation->reference_count++; - } - - iterator& operator=(const iterator& Iter) - { - if (Ptr_implementation != Iter.Ptr_implementation){ - if (Ptr_implementation != 0 && --(Ptr_implementation->reference_count)==0) { - delete Ptr_implementation; - } - Ptr_implementation = Iter.Ptr_implementation; - if (Ptr_implementation != 0) Ptr_implementation->reference_count++; - } - return *this; - } - - - const Point_with_transformed_distance& - operator* () const - { - return *(*Ptr_implementation); - } - - // -> operator - const Point_with_transformed_distance* - operator-> () const - { - return &*(*Ptr_implementation); - } - - // prefix operator - iterator& - operator++() - { - ++(*Ptr_implementation); - return *this; - } - - // postfix operator - Object_wrapper - operator++(int) - { - return (*Ptr_implementation)++; - } - - - bool - operator==(const iterator& It) const - { - if ( - ((Ptr_implementation == 0) || - Ptr_implementation->Item_PriorityQueue.empty()) && - ((It.Ptr_implementation == 0) || - It.Ptr_implementation->Item_PriorityQueue.empty()) - ) - return true; - // else - return (Ptr_implementation == It.Ptr_implementation); - } - - bool - operator!=(const iterator& It) const - { - return !(*this == It); - } - - std::ostream& - statistics (std::ostream& s) - { - Ptr_implementation->statistics(s); - return s; - } - - ~iterator() - { - if (Ptr_implementation != 0) { - Ptr_implementation->reference_count--; - if (Ptr_implementation->reference_count==0) { - delete Ptr_implementation; - Ptr_implementation = 0; - } - } - } - - - }; // class iterator - - //data members - const Tree& m_tree; - Query_item m_query; - Distance m_dist; - FT m_Eps; - bool m_search_nearest; - }; // class - - template - void swap (typename Orthogonal_incremental_neighbor_search::iterator& x, - typename Orthogonal_incremental_neighbor_search::iterator& y) - { - typename Orthogonal_incremental_neighbor_search::iterator::Iterator_implementation - *tmp = x.Ptr_implementation; - x.Ptr_implementation = y.Ptr_implementation; - y.Ptr_implementation = tmp; - } - -} // namespace CGAL - -#endif // CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH_H diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h index e9f455d7..253c89b4 100644 --- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -34,141 +34,141 @@ namespace persistence_diagram { * \ingroup bottleneck_distance */ class Graph_matching { -public: - /** \internal \brief Constructor constructing an empty matching. */ - explicit Graph_matching(Persistence_graph &g); - /** \internal \brief Copy operator. */ - Graph_matching& operator=(const Graph_matching& m); - /** \internal \brief Is the matching perfect ? */ - bool perfect() const; - /** \internal \brief Augments the matching with a maximal set of edge-disjoint shortest augmenting paths. */ - bool multi_augment(); - /** \internal \brief Sets the maximum length of the edges allowed to be added in the matching, 0 initially. */ - void set_r(double r); - -private: - Persistence_graph& g; - double r; - /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */ - std::vector v_to_u; - /** \internal \brief All the unmatched points in U. */ - std::list unmatched_in_u; - - /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ - Layered_neighbors_finder layering() const; - /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */ - bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth); - /** \internal \brief Update the matching with the simple augmenting path given as parameter. */ - void update(std::vector & path); + public: + /** \internal \brief Constructor constructing an empty matching. */ + explicit Graph_matching(Persistence_graph &g); + /** \internal \brief Copy operator. */ + Graph_matching& operator=(const Graph_matching& m); + /** \internal \brief Is the matching perfect ? */ + bool perfect() const; + /** \internal \brief Augments the matching with a maximal set of edge-disjoint shortest augmenting paths. */ + bool multi_augment(); + /** \internal \brief Sets the maximum length of the edges allowed to be added in the matching, 0 initially. */ + void set_r(double r); + + private: + Persistence_graph& g; + double r; + /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */ + std::vector v_to_u; + /** \internal \brief All the unmatched points in U. */ + std::list unmatched_in_u; + + /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ + Layered_neighbors_finder layering() const; + /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */ + bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth); + /** \internal \brief Update the matching with the simple augmenting path given as parameter. */ + void update(std::vector & path); }; inline Graph_matching::Graph_matching(Persistence_graph& g) : g(g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u() { - for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index) - unmatched_in_u.emplace_back(u_point_index); + for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index) + unmatched_in_u.emplace_back(u_point_index); } inline Graph_matching& Graph_matching::operator=(const Graph_matching& m) { - g = m.g; - r = m.r; - v_to_u = m.v_to_u; - unmatched_in_u = m.unmatched_in_u; - return *this; + g = m.g; + r = m.r; + v_to_u = m.v_to_u; + unmatched_in_u = m.unmatched_in_u; + return *this; } inline bool Graph_matching::perfect() const { - return unmatched_in_u.empty(); + return unmatched_in_u.empty(); } inline bool Graph_matching::multi_augment() { - if (perfect()) - return false; - Layered_neighbors_finder layered_nf(layering()); - int max_depth = layered_nf.vlayers_number()*2 - 1; - double rn = sqrt(g.size()); - // verification of a necessary criterion in order to shortcut if possible - if (max_depth <0 || (unmatched_in_u.size() > rn && max_depth >= rn)) - return false; - bool successful = false; - std::list tries(unmatched_in_u); - for (auto it = tries.cbegin(); it != tries.cend(); it++) - // 'augment' has side-effects which have to be always executed, don't change order - successful = augment(layered_nf, *it, max_depth) || successful; - return successful; + if (perfect()) + return false; + Layered_neighbors_finder layered_nf(layering()); + int max_depth = layered_nf.vlayers_number()*2 - 1; + double rn = sqrt(g.size()); + // verification of a necessary criterion in order to shortcut if possible + if (max_depth < 0 || (unmatched_in_u.size() > rn && max_depth >= rn)) + return false; + bool successful = false; + std::list tries(unmatched_in_u); + for (auto it = tries.cbegin(); it != tries.cend(); it++) + // 'augment' has side-effects which have to be always executed, don't change order + successful = augment(layered_nf, *it, max_depth) || successful; + return successful; } inline void Graph_matching::set_r(double r) { - this->r = r; + this->r = r; } inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) { - //V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices. - std::vector path; - path.emplace_back(u_start_index); - do { - if (static_cast(path.size()) > max_depth) { - path.pop_back(); - path.pop_back(); - } - if (path.empty()) - return false; - path.emplace_back(layered_nf.pull_near(path.back(), static_cast(path.size())/2)); - while (path.back() == null_point_index()) { - path.pop_back(); - path.pop_back(); - if (path.empty()) - return false; - path.pop_back(); - path.emplace_back(layered_nf.pull_near(path.back(), path.size() / 2)); - } - path.emplace_back(v_to_u.at(path.back())); - } while (path.back() != null_point_index()); - //if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex - path.pop_back(); - update(path); - return true; + // V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices. + std::vector path; + path.emplace_back(u_start_index); + do { + if (static_cast (path.size()) > max_depth) { + path.pop_back(); + path.pop_back(); + } + if (path.empty()) + return false; + path.emplace_back(layered_nf.pull_near(path.back(), static_cast (path.size()) / 2)); + while (path.back() == null_point_index()) { + path.pop_back(); + path.pop_back(); + if (path.empty()) + return false; + path.pop_back(); + path.emplace_back(layered_nf.pull_near(path.back(), path.size() / 2)); + } + path.emplace_back(v_to_u.at(path.back())); + } while (path.back() != null_point_index()); + // if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex + path.pop_back(); + update(path); + return true; } inline Layered_neighbors_finder Graph_matching::layering() const { - std::list u_vertices(unmatched_in_u); - std::list v_vertices; - Neighbors_finder nf(g, r); - for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index) - nf.add(v_point_index); - Layered_neighbors_finder layered_nf(g, r); - for(int layer = 0; !u_vertices.empty(); layer++) { - // one layer is one step in the BFS - for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) { - std::vector u_succ(nf.pull_all_near(*it1)); - for (auto it2 = u_succ.begin(); it2 != u_succ.end(); ++it2) { - layered_nf.add(*it2, layer); - v_vertices.emplace_back(*it2); - } - } - // When the above for finishes, we have progress of one half-step (from U to V) in the BFS - u_vertices.clear(); - bool end = false; - for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) - if (v_to_u.at(*it) == null_point_index()) - // we stop when a nearest exposed V vertex (from U exposed vertices) has been found - end = true; - else - u_vertices.emplace_back(v_to_u.at(*it)); - // When the above for finishes, we have progress of one half-step (from V to U) in the BFS - if (end) - return layered_nf; - v_vertices.clear(); + std::list u_vertices(unmatched_in_u); + std::list v_vertices; + Neighbors_finder nf(g, r); + for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index) + nf.add(v_point_index); + Layered_neighbors_finder layered_nf(g, r); + for (int layer = 0; !u_vertices.empty(); layer++) { + // one layer is one step in the BFS + for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) { + std::vector u_succ(nf.pull_all_near(*it1)); + for (auto it2 = u_succ.begin(); it2 != u_succ.end(); ++it2) { + layered_nf.add(*it2, layer); + v_vertices.emplace_back(*it2); + } } - return layered_nf; + // When the above for finishes, we have progress of one half-step (from U to V) in the BFS + u_vertices.clear(); + bool end = false; + for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) + if (v_to_u.at(*it) == null_point_index()) + // we stop when a nearest exposed V vertex (from U exposed vertices) has been found + end = true; + else + u_vertices.emplace_back(v_to_u.at(*it)); + // When the above for finishes, we have progress of one half-step (from V to U) in the BFS + if (end) + return layered_nf; + v_vertices.clear(); + } + return layered_nf; } inline void Graph_matching::update(std::vector& path) { - unmatched_in_u.remove(path.front()); - for (auto it = path.cbegin(); it != path.cend(); ++it) { - // Be careful, the iterator is incremented twice each time - int tmp = *it; - v_to_u[*(++it)] = tmp; - } + unmatched_in_u.remove(path.front()); + for (auto it = path.cbegin(); it != path.cend(); ++it) { + // Be careful, the iterator is incremented twice each time + int tmp = *it; + v_to_u[*(++it)] = tmp; + } } diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h index 70342d0c..0b2d26fe 100644 --- a/src/Bottleneck_distance/include/gudhi/Internal_point.h +++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -32,39 +32,60 @@ int null_point_index(); /** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ struct Internal_point { - double vec[2]; - int point_index; - Internal_point() {} - Internal_point(double x, double y, int p_i) { 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 point_index==p.point_index; - } - bool operator!=(const Internal_point& p) const { return !(*this == p); } + double vec[2]; + int point_index; + + Internal_point() { } + + Internal_point(double x, double y, int p_i) { + 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 point_index == p.point_index; + } + + bool operator!=(const Internal_point& p) const { + return !(*this == p); + } }; inline int null_point_index() { - return -1; + return -1; } struct Construct_coord_iterator { - typedef const double* result_type; - const double* operator()(const Internal_point& p) const - { return p.vec; } - const double* operator()(const Internal_point& p, int) const - { return p.vec+2; } + typedef const double* result_type; + + const double* operator()(const Internal_point& p) const { + return p.vec; + } + + const double* operator()(const Internal_point& p, int) const { + return p.vec + 2; + } }; } // namespace persistence_diagram } // namespace Gudhi - - - - #endif // INTERNAL_POINT_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h index 792925b7..96ece360 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -24,10 +24,9 @@ #define NEIGHBORS_FINDER_H_ // Inclusion order is important for CGAL patch -#include "CGAL/Kd_tree_node.h" -#include "CGAL/Kd_tree.h" -#include "CGAL/Orthogonal_k_neighbor_search.h" - +#include +#include +#include #include #include @@ -43,123 +42,126 @@ namespace persistence_diagram { /** \internal \brief data structure used to find any point (including projections) in V near to a query point from U * (which can be a projection). * - * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically removed. + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically + * removed. * * \ingroup bottleneck_distance */ class Neighbors_finder { - - typedef CGAL::Dimension_tag<2> D; - typedef CGAL::Search_traits Traits; - typedef CGAL::Weighted_Minkowski_distance Distance; - typedef CGAL::Orthogonal_k_neighbor_search K_neighbor_search; - typedef K_neighbor_search::Tree Kd_tree; - -public: - /** \internal \brief Constructor taking the near distance definition as parameter. */ - Neighbors_finder(const Persistence_graph& g, double r); - /** \internal \brief A point added will be possibly pulled. */ - void add(int v_point_index); - /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ - int pull_near(int u_point_index); - /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */ - std::vector pull_all_near(int u_point_index); - -private: - const Persistence_graph& g; - const double r; - Kd_tree kd_t; - std::unordered_set projections_f; + typedef CGAL::Dimension_tag<2> D; + typedef CGAL::Search_traits Traits; + typedef CGAL::Weighted_Minkowski_distance Distance; + typedef CGAL::Orthogonal_k_neighbor_search K_neighbor_search; + typedef K_neighbor_search::Tree Kd_tree; + + public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if + * there isn't such a point. */ + int pull_near(int u_point_index); + /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */ + std::vector pull_all_near(int u_point_index); + + private: + const Persistence_graph& g; + const double r; + Kd_tree kd_t; + std::unordered_set projections_f; }; /** \internal \brief data structure used to find any point (including projections) in V near to a query point from U * (which can be a projection) in a layered graph layer given as parmeter. * - * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically removed. + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically + * removed. * * \ingroup bottleneck_distance */ class Layered_neighbors_finder { -public: - /** \internal \brief Constructor taking the near distance definition as parameter. */ - Layered_neighbors_finder(const Persistence_graph& g, double r); - /** \internal \brief A point added will be possibly pulled. */ - void add(int v_point_index, int vlayer); - /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if there isn't such a point. */ - int pull_near(int u_point_index, int vlayer); - /** \internal \brief Returns the number of layers. */ - int vlayers_number() const; - -private: - const Persistence_graph& g; - const double r; - std::vector> neighbors_finder; + public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Layered_neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index, int vlayer); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if + * there isn't such a point. */ + int pull_near(int u_point_index, int vlayer); + /** \internal \brief Returns the number of layers. */ + int vlayers_number() const; + + private: + const Persistence_graph& g; + const double r; + std::vector> neighbors_finder; }; inline Neighbors_finder::Neighbors_finder(const Persistence_graph& g, double r) : g(g), r(r), kd_t(), projections_f() { } inline void Neighbors_finder::add(int v_point_index) { - if (g.on_the_v_diagonal(v_point_index)) - projections_f.emplace(v_point_index); - else - kd_t.insert(g.get_v_point(v_point_index)); + if (g.on_the_v_diagonal(v_point_index)) + projections_f.emplace(v_point_index); + else + kd_t.insert(g.get_v_point(v_point_index)); } inline int Neighbors_finder::pull_near(int u_point_index) { - int tmp; - int c = g.corresponding_point_in_v(u_point_index); - if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()){ - //Any pair of projection is at distance 0 - tmp = *projections_f.cbegin(); - projections_f.erase(tmp); - } - else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)){ - //Is the query point near to its projection ? - tmp = c; - projections_f.erase(tmp); - } - else{ - //Is the query point near to a V point in the plane ? - Internal_point u_point = g.get_u_point(u_point_index); - std::array w = { {1., 1.} }; - K_neighbor_search search(kd_t, u_point, 1, 0., true, Distance(0, 2, w.begin(), w.end())); - auto it = search.begin(); - if(it==search.end() || g.distance(u_point_index, it->first.point_index) > r) - return null_point_index(); - tmp = it->first.point_index; - kd_t.remove(g.get_v_point(tmp)); - } - return tmp; + int tmp; + int c = g.corresponding_point_in_v(u_point_index); + if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()) { + // Any pair of projection is at distance 0 + tmp = *projections_f.cbegin(); + projections_f.erase(tmp); + } else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)) { + // Is the query point near to its projection ? + tmp = c; + projections_f.erase(tmp); + } else { + // Is the query point near to a V point in the plane ? + Internal_point u_point = g.get_u_point(u_point_index); + std::array w = { + {1., 1.} + }; + K_neighbor_search search(kd_t, u_point, 1, 0., true, Distance(0, 2, w.begin(), w.end())); + auto it = search.begin(); + if (it == search.end() || g.distance(u_point_index, it->first.point_index) > r) + return null_point_index(); + tmp = it->first.point_index; + kd_t.remove(g.get_v_point(tmp)); + } + return tmp; } inline std::vector Neighbors_finder::pull_all_near(int u_point_index) { - std::vector all_pull; - int last_pull = pull_near(u_point_index); - while (last_pull != null_point_index()) { - all_pull.emplace_back(last_pull); - last_pull = pull_near(u_point_index); - } - return all_pull; + std::vector all_pull; + int last_pull = pull_near(u_point_index); + while (last_pull != null_point_index()) { + all_pull.emplace_back(last_pull); + last_pull = pull_near(u_point_index); + } + return all_pull; } inline Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_graph& g, double r) : g(g), r(r), neighbors_finder() { } inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) { - for (int l = neighbors_finder.size(); l <= vlayer; l++) - neighbors_finder.emplace_back(std::unique_ptr(new Neighbors_finder(g, r))); - neighbors_finder.at(vlayer)->add(v_point_index); + for (int l = neighbors_finder.size(); l <= vlayer; l++) + neighbors_finder.emplace_back(std::unique_ptr(new Neighbors_finder(g, r))); + neighbors_finder.at(vlayer)->add(v_point_index); } inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) { - if (static_cast (neighbors_finder.size()) <= vlayer) - return null_point_index(); - return neighbors_finder.at(vlayer)->pull_near(u_point_index); + if (static_cast (neighbors_finder.size()) <= vlayer) + return null_point_index(); + return neighbors_finder.at(vlayer)->pull_near(u_point_index); } inline int Layered_neighbors_finder::vlayers_number() const { - return static_cast(neighbors_finder.size()); + return static_cast (neighbors_finder.size()); } } // namespace persistence_diagram diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h index 45a4d586..3a4a5fec 100644 --- a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h +++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -31,147 +31,144 @@ namespace Gudhi { namespace persistence_diagram { - /** \internal \brief Structure representing an euclidean bipartite graph containing * the points from the two persistence diagrams (including the projections). * * \ingroup bottleneck_distance */ class Persistence_graph { -public: - /** \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 ? */ - bool on_the_u_diagonal(int u_point_index) const; - /** \internal \brief Is the given point from V the projection of a point in U ? */ - bool on_the_v_diagonal(int v_point_index) const; - /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */ - int corresponding_point_in_u(int v_point_index) const; - /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */ - int corresponding_point_in_v(int u_point_index) const; - /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */ - 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 ? */ - double bottleneck_alive() const; - /** \internal \brief Returns the O(n^2) sorted distances between the points. */ - std::vector sorted_distances() const; - /** \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; - /** \internal \brief Returns the corresponding internal point */ - Internal_point get_v_point(int v_point_index) const; - -private: - std::vector u; - std::vector v; - double b_alive; + public: + /** \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 ? */ + bool on_the_u_diagonal(int u_point_index) const; + /** \internal \brief Is the given point from V the projection of a point in U ? */ + bool on_the_v_diagonal(int v_point_index) const; + /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */ + int corresponding_point_in_u(int v_point_index) const; + /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */ + int corresponding_point_in_v(int u_point_index) const; + /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */ + 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 ? */ + double bottleneck_alive() const; + /** \internal \brief Returns the O(n^2) sorted distances between the points. */ + std::vector sorted_distances() const; + /** \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; + /** \internal \brief Returns the corresponding internal point */ + Internal_point get_v_point(int v_point_index) const; + + private: + std::vector u; + std::vector v; + double b_alive; }; template Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e) - : u(), v(), b_alive(0.) -{ - std::vector u_alive; - std::vector v_alive; - for (auto it = std::begin(diag1); it != std::end(diag1); ++it){ - if(std::get<1>(*it) == std::numeric_limits::infinity()) - 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::infinity()) - 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())); - } - if (u.size() < v.size()) - swap(u, v); - std::sort(u_alive.begin(), u_alive.end()); - 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) - b_alive = std::max(b_alive, std::fabs(*it_u - *it_v)); + : u(), v(), b_alive(0.) { + std::vector u_alive; + std::vector v_alive; + for (auto it = std::begin(diag1); it != std::end(diag1); ++it) { + if (std::get<1>(*it) == std::numeric_limits::infinity()) + 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::infinity()) + 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())); + } + if (u.size() < v.size()) + swap(u, v); + std::sort(u_alive.begin(), u_alive.end()); + 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) + b_alive = std::max(b_alive, std::fabs(*it_u - *it_v)); } inline bool Persistence_graph::on_the_u_diagonal(int u_point_index) const { - return u_point_index >= static_cast (u.size()); + return u_point_index >= static_cast (u.size()); } inline bool Persistence_graph::on_the_v_diagonal(int v_point_index) const { - return v_point_index >= static_cast (v.size()); + return v_point_index >= static_cast (v.size()); } inline int Persistence_graph::corresponding_point_in_u(int v_point_index) const { - return on_the_v_diagonal(v_point_index) ? - v_point_index - static_cast (v.size()) : v_point_index + static_cast (u.size()); + return on_the_v_diagonal(v_point_index) ? + v_point_index - static_cast (v.size()) : v_point_index + static_cast (u.size()); } inline int Persistence_graph::corresponding_point_in_v(int u_point_index) const { - return on_the_u_diagonal(u_point_index) ? - u_point_index - static_cast (u.size()) : u_point_index + static_cast (v.size()); + return on_the_u_diagonal(u_point_index) ? + u_point_index - static_cast (u.size()) : u_point_index + static_cast (v.size()); } inline double Persistence_graph::distance(int u_point_index, int v_point_index) const { - if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(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); - return std::max(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y())); + if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(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); + return std::max(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y())); } inline int Persistence_graph::size() const { - return static_cast (u.size() + v.size()); + return static_cast (u.size() + v.size()); } -inline double Persistence_graph::bottleneck_alive() const{ - return b_alive; +inline double Persistence_graph::bottleneck_alive() const { + return b_alive; } inline std::vector Persistence_graph::sorted_distances() const { - std::vector distances; - 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) - distances.push_back(distance(u_point_index, v_point_index)); - } - std::sort(distances.begin(), distances.end()); - return distances; + std::vector distances; + 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) + 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 { - 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.; - return Internal_point(m,m,u_point_index); + 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.; + return Internal_point(m, m, u_point_index); } 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.; - return Internal_point(m,m,v_point_index); + 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.; + return Internal_point(m, m, v_point_index); } 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()); - for(auto it = v.cbegin(); it != v.cend(); it++) - max = std::max(max, it->y()); - return max; + double max = 0.; + for (auto it = u.cbegin(); it != u.cend(); it++) + max = std::max(max, it->y()); + for (auto it = v.cbegin(); it != v.cend(); it++) + max = std::max(max, it->y()); + return max; } - } // namespace persistence_diagram } // namespace Gudhi diff --git a/src/Bottleneck_distance/test/CMakeLists.txt b/src/Bottleneck_distance/test/CMakeLists.txt index 13213075..a6979d3c 100644 --- a/src/Bottleneck_distance/test/CMakeLists.txt +++ b/src/Bottleneck_distance/test/CMakeLists.txt @@ -17,7 +17,6 @@ if(CGAL_FOUND) if (NOT CGAL_VERSION VERSION_LESS 4.8.0) if (EIGEN3_FOUND) add_executable ( bottleneckUT bottleneck_unit_test.cpp ) - add_executable ( bottleneck_chrono bottleneck_chrono.cpp ) target_link_libraries(bottleneckUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) # Unitary tests diff --git a/src/Bottleneck_distance/test/bottleneck_chrono.cpp b/src/Bottleneck_distance/test/bottleneck_chrono.cpp deleted file mode 100644 index a30d42b5..00000000 --- a/src/Bottleneck_distance/test/bottleneck_chrono.cpp +++ /dev/null @@ -1,62 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author: Francois Godi - * - * Copyright (C) 2015 INRIA (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#include -#include -#include -#include - -using namespace Gudhi::persistence_diagram; - - -double upper_bound = 400.; // any real >0 - -int main(){ - std::ofstream objetfichier; - objetfichier.open("results.csv", std::ios::out); - - for(int n = 1000; n<=10000; n+=1000){ - std::uniform_real_distribution unif1(0.,upper_bound); - std::uniform_real_distribution unif2(upper_bound/1000.,upper_bound/100.); - std::default_random_engine re; - std::vector< std::pair > v1, v2; - for (int i = 0; i < n; 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); - } - std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now(); - double b = bottleneck_distance(v1, v2); - 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)); - objetfichier << n << ";" << duration.count() << ";" << b << std::endl; - } - objetfichier.close(); -} diff --git a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp index fba1d369..e39613b3 100644 --- a/src/Bottleneck_distance/test/bottleneck_unit_test.cpp +++ b/src/Bottleneck_distance/test/bottleneck_unit_test.cpp @@ -4,7 +4,7 @@ * * Author: Francois Godi * - * Copyright (C) 2015 INRIA (France) + * Copyright (C) 2015 INRIA * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -30,138 +30,138 @@ using namespace Gudhi::persistence_diagram; -int n1 = 81; // a natural number >0 -int n2 = 180; // a natural number >0 -double upper_bound = 406.43; // any real >0 +int n1 = 81; // a natural number >0 +int n2 = 180; // a natural number >0 +double upper_bound = 406.43; // any real >0 -std::uniform_real_distribution unif(0.,upper_bound); +std::uniform_real_distribution unif(0., upper_bound); std::default_random_engine re; std::vector< std::pair > v1, v2; -BOOST_AUTO_TEST_CASE(persistence_graph){ - // Random construction - for (int i = 0; i < n1; i++) { - double a = unif(re); - double b = unif(re); - v1.emplace_back(std::min(a,b), std::max(a,b)); - } - for (int i = 0; i < n2; i++) { - double a = unif(re); - double b = unif(re); - v2.emplace_back(std::min(a,b), std::max(a,b)); - } - Persistence_graph g(v1, v2, 0.); - std::vector d(g.sorted_distances()); - // - BOOST_CHECK(!g.on_the_u_diagonal(n1-1)); - BOOST_CHECK(!g.on_the_u_diagonal(n1)); - BOOST_CHECK(!g.on_the_u_diagonal(n2-1)); - BOOST_CHECK(g.on_the_u_diagonal(n2)); - BOOST_CHECK(!g.on_the_v_diagonal(n1-1)); - BOOST_CHECK(g.on_the_v_diagonal(n1)); - BOOST_CHECK(g.on_the_v_diagonal(n2-1)); - BOOST_CHECK(g.on_the_v_diagonal(n2)); - // - BOOST_CHECK(g.corresponding_point_in_u(0)==n2); - BOOST_CHECK(g.corresponding_point_in_u(n1)==0); - BOOST_CHECK(g.corresponding_point_in_v(0)==n1); - BOOST_CHECK(g.corresponding_point_in_v(n2)==0); - // - BOOST_CHECK(g.size()==(n1+n2)); - // - 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(persistence_graph) { + // Random construction + for (int i = 0; i < n1; i++) { + double a = unif(re); + double b = unif(re); + v1.emplace_back(std::min(a, b), std::max(a, b)); + } + for (int i = 0; i < n2; i++) { + double a = unif(re); + double b = unif(re); + v2.emplace_back(std::min(a, b), std::max(a, b)); + } + Persistence_graph g(v1, v2, 0.); + std::vector d(g.sorted_distances()); + // + BOOST_CHECK(!g.on_the_u_diagonal(n1 - 1)); + BOOST_CHECK(!g.on_the_u_diagonal(n1)); + BOOST_CHECK(!g.on_the_u_diagonal(n2 - 1)); + BOOST_CHECK(g.on_the_u_diagonal(n2)); + BOOST_CHECK(!g.on_the_v_diagonal(n1 - 1)); + BOOST_CHECK(g.on_the_v_diagonal(n1)); + BOOST_CHECK(g.on_the_v_diagonal(n2 - 1)); + BOOST_CHECK(g.on_the_v_diagonal(n2)); + // + BOOST_CHECK(g.corresponding_point_in_u(0) == n2); + BOOST_CHECK(g.corresponding_point_in_u(n1) == 0); + BOOST_CHECK(g.corresponding_point_in_v(0) == n1); + BOOST_CHECK(g.corresponding_point_in_v(n2) == 0); + // + BOOST_CHECK(g.size() == (n1 + n2)); + // + 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) { - Persistence_graph g(v1, v2, 0.); - Neighbors_finder nf(g, 1.); - for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2) - nf.add(v_point_index); - // - int v_point_index_1 = nf.pull_near(n2/2); - BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); - std::vector l = nf.pull_all_near(n2/2); - bool v = true; - for(auto it = l.cbegin(); it != l.cend(); ++it) - v = v && (g.distance(n2/2,*it)>1.); - BOOST_CHECK(v); - int v_point_index_2 = nf.pull_near(n2/2); - BOOST_CHECK(v_point_index_2 == -1); + Persistence_graph g(v1, v2, 0.); + Neighbors_finder nf(g, 1.); + for (int v_point_index = 1; v_point_index < ((n2 + n1)*9 / 10); v_point_index += 2) + nf.add(v_point_index); + // + int v_point_index_1 = nf.pull_near(n2 / 2); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + std::vector l = nf.pull_all_near(n2 / 2); + bool v = true; + for (auto it = l.cbegin(); it != l.cend(); ++it) + v = v && (g.distance(n2 / 2, *it) > 1.); + BOOST_CHECK(v); + int v_point_index_2 = nf.pull_near(n2 / 2); + BOOST_CHECK(v_point_index_2 == -1); } BOOST_AUTO_TEST_CASE(layered_neighbors_finder) { - Persistence_graph g(v1, v2, 0.); - Layered_neighbors_finder lnf(g, 1.); - for(int v_point_index=1; v_point_index<((n2+n1)*9/10); v_point_index+=2) - lnf.add(v_point_index, v_point_index % 7); - // - int v_point_index_1 = lnf.pull_near(n2/2,6); - BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); - int v_point_index_2 = lnf.pull_near(n2/2,6); - BOOST_CHECK(v_point_index_2 == -1); - v_point_index_1 = lnf.pull_near(n2/2,0); - BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2/2,v_point_index_1)<=1.)); - v_point_index_2 = lnf.pull_near(n2/2,0); - BOOST_CHECK(v_point_index_2 == -1); + Persistence_graph g(v1, v2, 0.); + Layered_neighbors_finder lnf(g, 1.); + for (int v_point_index = 1; v_point_index < ((n2 + n1)*9 / 10); v_point_index += 2) + lnf.add(v_point_index, v_point_index % 7); + // + int v_point_index_1 = lnf.pull_near(n2 / 2, 6); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + int v_point_index_2 = lnf.pull_near(n2 / 2, 6); + BOOST_CHECK(v_point_index_2 == -1); + v_point_index_1 = lnf.pull_near(n2 / 2, 0); + BOOST_CHECK((v_point_index_1 == -1) || (g.distance(n2 / 2, v_point_index_1) <= 1.)); + v_point_index_2 = lnf.pull_near(n2 / 2, 0); + BOOST_CHECK(v_point_index_2 == -1); } BOOST_AUTO_TEST_CASE(graph_matching) { - Persistence_graph g(v1, v2, 0.); - Graph_matching m1(g); - m1.set_r(0.); - int e = 0; - while (m1.multi_augment()) - ++e; - BOOST_CHECK(e > 0); - BOOST_CHECK(e <= 2*sqrt(2*(n1+n2))); - Graph_matching m2 = m1; - BOOST_CHECK(!m2.multi_augment()); - m2.set_r(upper_bound); - e = 0; - while (m2.multi_augment()) - ++e; - BOOST_CHECK(e <= 2*sqrt(2*(n1+n2))); - BOOST_CHECK(m2.perfect()); - BOOST_CHECK(!m1.perfect()); + Persistence_graph g(v1, v2, 0.); + Graph_matching m1(g); + m1.set_r(0.); + int e = 0; + while (m1.multi_augment()) + ++e; + BOOST_CHECK(e > 0); + BOOST_CHECK(e <= 2 * sqrt(2 * (n1 + n2))); + Graph_matching m2 = m1; + BOOST_CHECK(!m2.multi_augment()); + m2.set_r(upper_bound); + e = 0; + while (m2.multi_augment()) + ++e; + BOOST_CHECK(e <= 2 * sqrt(2 * (n1 + n2))); + BOOST_CHECK(m2.perfect()); + BOOST_CHECK(!m1.perfect()); } -BOOST_AUTO_TEST_CASE(global){ - std::uniform_real_distribution unif1(0.,upper_bound); - 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++) { - 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, 0.) <= upper_bound/100.); - BOOST_CHECK(bottleneck_distance(v1, v2, upper_bound/10000.) <= upper_bound/100. + upper_bound/10000.); - BOOST_CHECK(std::abs(bottleneck_distance(v1, v2, 0.) - bottleneck_distance(v1, v2, upper_bound/10000.)) <= upper_bound/10000.); +BOOST_AUTO_TEST_CASE(global) { + std::uniform_real_distribution unif1(0., upper_bound); + 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++) { + 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, 0.) <= upper_bound / 100.); + BOOST_CHECK(bottleneck_distance(v1, v2, upper_bound / 10000.) <= upper_bound / 100. + upper_bound / 10000.); + BOOST_CHECK(std::abs(bottleneck_distance(v1, v2, 0.) - bottleneck_distance(v1, v2, upper_bound / 10000.)) <= upper_bound / 10000.); } diff --git a/src/Tangential_complex/benchmark/CMakeLists.txt b/src/Tangential_complex/benchmark/CMakeLists.txt index 56dd8128..788c2b4d 100644 --- a/src/Tangential_complex/benchmark/CMakeLists.txt +++ b/src/Tangential_complex/benchmark/CMakeLists.txt @@ -1,15 +1,6 @@ cmake_minimum_required(VERSION 2.6) project(Tangential_complex_benchmark) -if (GCOVR_PATH) - # for gcovr to make coverage reports - Corbera Jenkins plugin - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") -endif() -if (GPROF_PATH) - # for gprof to make coverage reports - Jenkins - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") -endif() - # need CGAL 4.8 if(CGAL_FOUND) if (NOT CGAL_VERSION VERSION_LESS 4.8.0) diff --git a/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt b/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt new file mode 100644 index 00000000..a588d113 --- /dev/null +++ b/src/common/include/gudhi_patches/Bottleneck_distance_CGAL_patches.txt @@ -0,0 +1,3 @@ +CGAL/Kd_tree.h +CGAL/Kd_tree_node.h +CGAL/Orthogonal_incremental_neighbor_search.h diff --git a/src/common/include/gudhi_patches/CGAL/Kd_tree.h b/src/common/include/gudhi_patches/CGAL/Kd_tree.h new file mode 100644 index 00000000..f085b0da --- /dev/null +++ b/src/common/include/gudhi_patches/CGAL/Kd_tree.h @@ -0,0 +1,582 @@ +// Copyright (c) 2002,2011,2014 Utrecht University (The Netherlands), Max-Planck-Institute Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Hans Tangelder (), +// : Waqar Khan + +#ifndef CGAL_KD_TREE_H +#define CGAL_KD_TREE_H + +#include "Kd_tree_node.h" + +#include +#include +#include + +#include +#include +#include + + +#include +#include +#include + +#ifdef CGAL_HAS_THREADS +#include +#endif + +namespace CGAL { + +//template , class UseExtendedNode = Tag_true > +template , class UseExtendedNode = Tag_true > +class Kd_tree { + +public: + typedef SearchTraits Traits; + typedef Splitter_ Splitter; + typedef typename SearchTraits::Point_d Point_d; + typedef typename Splitter::Container Point_container; + + typedef typename SearchTraits::FT FT; + typedef Kd_tree_node Node; + typedef Kd_tree_leaf_node Leaf_node; + typedef Kd_tree_internal_node Internal_node; + typedef Kd_tree Tree; + typedef Kd_tree Self; + + typedef Node* Node_handle; + typedef const Node* Node_const_handle; + typedef Leaf_node* Leaf_node_handle; + typedef const Leaf_node* Leaf_node_const_handle; + typedef Internal_node* Internal_node_handle; + typedef const Internal_node* Internal_node_const_handle; + typedef typename std::vector::const_iterator Point_d_iterator; + typedef typename std::vector::const_iterator Point_d_const_iterator; + typedef typename Splitter::Separator Separator; + typedef typename std::vector::const_iterator iterator; + typedef typename std::vector::const_iterator const_iterator; + + typedef typename std::vector::size_type size_type; + + typedef typename internal::Get_dimension_tag::Dimension D; + +private: + SearchTraits traits_; + Splitter split; + + + // wokaround for https://svn.boost.org/trac/boost/ticket/9332 +#if (_MSC_VER == 1800) && (BOOST_VERSION == 105500) + std::deque internal_nodes; + std::deque leaf_nodes; +#else + boost::container::deque internal_nodes; + boost::container::deque leaf_nodes; +#endif + + Node_handle tree_root; + + Kd_tree_rectangle* bbox; + std::vector pts; + + // Instead of storing the points in arrays in the Kd_tree_node + // we put all the data in a vector in the Kd_tree. + // and we only store an iterator range in the Kd_tree_node. + // + std::vector data; + + + #ifdef CGAL_HAS_THREADS + mutable CGAL_MUTEX building_mutex;//mutex used to protect const calls inducing build() + #endif + bool built_; + bool removed_; + + // protected copy constructor + Kd_tree(const Tree& tree) + : traits_(tree.traits_),built_(tree.built_) + {}; + + + // Instead of the recursive construction of the tree in the class Kd_tree_node + // we do this in the tree class. The advantage is that we then can optimize + // the allocation of the nodes. + + // The leaf node + Node_handle + create_leaf_node(Point_container& c) + { + Leaf_node node(true , static_cast(c.size())); + std::ptrdiff_t tmp = c.begin() - data.begin(); + node.data = pts.begin() + tmp; + + leaf_nodes.push_back(node); + Leaf_node_handle nh = &leaf_nodes.back(); + + + return nh; + } + + + // The internal node + + Node_handle + create_internal_node(Point_container& c, const Tag_true&) + { + return create_internal_node_use_extension(c); + } + + Node_handle + create_internal_node(Point_container& c, const Tag_false&) + { + return create_internal_node(c); + } + + + + // TODO: Similiar to the leaf_init function above, a part of the code should be + // moved to a the class Kd_tree_node. + // It is not proper yet, but the goal was to see if there is + // a potential performance gain through the Compact_container + Node_handle + create_internal_node_use_extension(Point_container& c) + { + Internal_node node(false); + internal_nodes.push_back(node); + Internal_node_handle nh = &internal_nodes.back(); + + Separator sep; + Point_container c_low(c.dimension(),traits_); + split(sep, c, c_low); + nh->set_separator(sep); + + int cd = nh->cutting_dimension(); + if(!c_low.empty()){ + nh->lower_low_val = c_low.tight_bounding_box().min_coord(cd); + nh->lower_high_val = c_low.tight_bounding_box().max_coord(cd); + } + else{ + nh->lower_low_val = nh->cutting_value(); + nh->lower_high_val = nh->cutting_value(); + } + if(!c.empty()){ + nh->upper_low_val = c.tight_bounding_box().min_coord(cd); + nh->upper_high_val = c.tight_bounding_box().max_coord(cd); + } + else{ + nh->upper_low_val = nh->cutting_value(); + nh->upper_high_val = nh->cutting_value(); + } + + CGAL_assertion(nh->cutting_value() >= nh->lower_low_val); + CGAL_assertion(nh->cutting_value() <= nh->upper_high_val); + + if (c_low.size() > split.bucket_size()){ + nh->lower_ch = create_internal_node_use_extension(c_low); + }else{ + nh->lower_ch = create_leaf_node(c_low); + } + if (c.size() > split.bucket_size()){ + nh->upper_ch = create_internal_node_use_extension(c); + }else{ + nh->upper_ch = create_leaf_node(c); + } + + + + + return nh; + } + + + // Note also that I duplicated the code to get rid if the if's for + // the boolean use_extension which was constant over the construction + Node_handle + create_internal_node(Point_container& c) + { + Internal_node node(false); + internal_nodes.push_back(node); + Internal_node_handle nh = &internal_nodes.back(); + Separator sep; + + Point_container c_low(c.dimension(),traits_); + split(sep, c, c_low); + nh->set_separator(sep); + + if (c_low.size() > split.bucket_size()){ + nh->lower_ch = create_internal_node(c_low); + }else{ + nh->lower_ch = create_leaf_node(c_low); + } + if (c.size() > split.bucket_size()){ + nh->upper_ch = create_internal_node(c); + }else{ + nh->upper_ch = create_leaf_node(c); + } + + + + return nh; + } + + + +public: + + Kd_tree(Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) + : traits_(traits),split(s), built_(false), removed_(false) + {} + + template + Kd_tree(InputIterator first, InputIterator beyond, + Splitter s = Splitter(),const SearchTraits traits=SearchTraits()) + : traits_(traits),split(s), built_(false), removed_(false) + { + pts.insert(pts.end(), first, beyond); + } + + bool empty() const { + return pts.empty(); + } + + void + build() + { + // This function is not ready to be called when a tree already exists, one + // must call invalidate_built() first. + CGAL_assertion(!is_built()); + CGAL_assertion(!removed_); + const Point_d& p = *pts.begin(); + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits_.construct_cartesian_const_iterator_d_object(); + int dim = static_cast(std::distance(ccci(p), ccci(p,0))); + + data.reserve(pts.size()); + for(unsigned int i = 0; i < pts.size(); i++){ + data.push_back(&pts[i]); + } + Point_container c(dim, data.begin(), data.end(),traits_); + bbox = new Kd_tree_rectangle(c.bounding_box()); + if (c.size() <= split.bucket_size()){ + tree_root = create_leaf_node(c); + }else { + tree_root = create_internal_node(c, UseExtendedNode()); + } + + //Reorder vector for spatial locality + std::vector ptstmp; + ptstmp.resize(pts.size()); + for (std::size_t i = 0; i < pts.size(); ++i){ + ptstmp[i] = *data[i]; + } + for(std::size_t i = 0; i < leaf_nodes.size(); ++i){ + std::ptrdiff_t tmp = leaf_nodes[i].begin() - pts.begin(); + leaf_nodes[i].data = ptstmp.begin() + tmp; + } + pts.swap(ptstmp); + + data.clear(); + + built_ = true; + } + +private: + //any call to this function is for the moment not threadsafe + void const_build() const { + #ifdef CGAL_HAS_THREADS + //this ensure that build() will be called once + CGAL_SCOPED_LOCK(building_mutex); + if(!is_built()) + #endif + const_cast(this)->build(); //THIS IS NOT THREADSAFE + } +public: + + bool is_built() const + { + return built_; + } + + void invalidate_built() + { + if(removed_){ + // Walk the tree to collect the remaining points. + // Writing directly to pts would likely work, but better be safe. + std::vector ptstmp; + //ptstmp.resize(root()->num_items()); + root()->tree_items(std::back_inserter(ptstmp)); + pts.swap(ptstmp); + removed_=false; + CGAL_assertion(is_built()); // the rest of the cleanup must happen + } + if(is_built()){ + internal_nodes.clear(); + leaf_nodes.clear(); + data.clear(); + delete bbox; + built_ = false; + } + } + + void clear() + { + invalidate_built(); + pts.clear(); + removed_ = false; + } + + void + insert(const Point_d& p) + { + invalidate_built(); + pts.push_back(p); + } + + template + void + insert(InputIterator first, InputIterator beyond) + { + invalidate_built(); + pts.insert(pts.end(),first, beyond); + } + +private: + struct Equal_by_coordinates { + SearchTraits const* traits; + Point_d const* pp; + bool operator()(Point_d const&q) const { + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits->construct_cartesian_const_iterator_d_object(); + return std::equal(ccci(*pp), ccci(*pp,0), ccci(q)); + } + }; + Equal_by_coordinates equal_by_coordinates(Point_d const&p){ + Equal_by_coordinates ret = { &traits(), &p }; + return ret; + } + +public: + void + remove(const Point_d& p) + { + remove(p, equal_by_coordinates(p)); + } + + template + void + remove(const Point_d& p, Equal const& equal_to_p) + { +#if 0 + // This code could have quadratic runtime. + if (!is_built()) { + std::vector::iterator pi = std::find(pts.begin(), pts.end(), p); + // Precondition: the point must be there. + CGAL_assertion (pi != pts.end()); + pts.erase(pi); + return; + } +#endif + bool success = remove_(p, 0, false, 0, false, root(), equal_to_p); + CGAL_assertion(success); + + // Do not set the flag is the tree has been cleared. + if(is_built()) + removed_ |= success; + } +private: + template + bool remove_(const Point_d& p, + Internal_node_handle grandparent, bool parent_islower, + Internal_node_handle parent, bool islower, + Node_handle node, Equal const& equal_to_p) { + // Recurse to locate the point + if (!node->is_leaf()) { + Internal_node_handle newparent = static_cast(node); + // FIXME: This should be if(xcutting_dimension()] <= newparent->cutting_value()) { + if (remove_(p, parent, islower, newparent, true, newparent->lower(), equal_to_p)) + return true; + } + //if (traits().construct_cartesian_const_iterator_d_object()(p)[newparent->cutting_dimension()] >= newparent->cutting_value()) + return remove_(p, parent, islower, newparent, false, newparent->upper(), equal_to_p); + + CGAL_assertion(false); // Point was not found + } + + // Actual removal + Leaf_node_handle lnode = static_cast(node); + if (lnode->size() > 1) { + iterator pi = std::find_if(lnode->begin(), lnode->end(), equal_to_p); + // FIXME: we should ensure this never happens + if (pi == lnode->end()) return false; + iterator lasti = lnode->end() - 1; + if (pi != lasti) { + // Hack to get a non-const iterator + std::iter_swap(pts.begin()+(pi-pts.begin()), pts.begin()+(lasti-pts.begin())); + } + lnode->drop_last_point(); + } else if (!equal_to_p(*lnode->begin())) { + // FIXME: we should ensure this never happens + return false; + } else if (grandparent) { + Node_handle brother = islower ? parent->upper() : parent->lower(); + if (parent_islower) + grandparent->set_lower(brother); + else + grandparent->set_upper(brother); + } else if (parent) { + tree_root = islower ? parent->upper() : parent->lower(); + } else { + clear(); + } + return true; + } + +public: + //For efficiency; reserve the size of the points vectors in advance (if the number of points is already known). + void reserve(size_t size) + { + pts.reserve(size); + } + + //Get the capacity of the underlying points vector. + size_t capacity() + { + return pts.capacity(); + } + + + template + OutputIterator + search(OutputIterator it, const FuzzyQueryItem& q) const + { + if(! pts.empty()){ + + if(! is_built()){ + const_build(); + } + Kd_tree_rectangle b(*bbox); + return tree_root->search(it,q,b); + } + return it; + } + + + template + boost::optional + search_any_point(const FuzzyQueryItem& q) const + { + if(! pts.empty()){ + + if(! is_built()){ + const_build(); + } + Kd_tree_rectangle b(*bbox); + return tree_root->search_any_point(q,b); + } + return boost::none; + } + + + ~Kd_tree() { + if(is_built()){ + delete bbox; + } + } + + + const SearchTraits& + traits() const + { + return traits_; + } + + Node_const_handle + root() const + { + if(! is_built()){ + const_build(); + } + return tree_root; + } + + Node_handle + root() + { + if(! is_built()){ + build(); + } + return tree_root; + } + + void + print() const + { + if(! is_built()){ + const_build(); + } + root()->print(); + } + + const Kd_tree_rectangle& + bounding_box() const + { + if(! is_built()){ + const_build(); + } + return *bbox; + } + + const_iterator + begin() const + { + return pts.begin(); + } + + const_iterator + end() const + { + return pts.end(); + } + + size_type + size() const + { + return pts.size(); + } + + // Print statistics of the tree. + std::ostream& + statistics(std::ostream& s) const + { + if(! is_built()){ + const_build(); + } + s << "Tree statistics:" << std::endl; + s << "Number of items stored: " + << root()->num_items() << std::endl; + s << "Number of nodes: " + << root()->num_nodes() << std::endl; + s << " Tree depth: " << root()->depth() << std::endl; + return s; + } + + +}; + +} // namespace CGAL + +#endif // CGAL_KD_TREE_H diff --git a/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h b/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h new file mode 100644 index 00000000..909ee260 --- /dev/null +++ b/src/common/include/gudhi_patches/CGAL/Kd_tree_node.h @@ -0,0 +1,586 @@ +// Copyright (c) 2002,2011 Utrecht University (The Netherlands). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Authors : Hans Tangelder () + +#ifndef CGAL_KD_TREE_NODE_H +#define CGAL_KD_TREE_NODE_H + +#include "CGAL/Splitters.h" + +#include +#include + +namespace CGAL { + + template + class Kd_tree; + + template < class TreeTraits, class Splitter, class UseExtendedNode > + class Kd_tree_node { + + friend class Kd_tree; + + typedef typename Kd_tree::Node_handle Node_handle; + typedef typename Kd_tree::Node_const_handle Node_const_handle; + typedef typename Kd_tree::Internal_node_handle Internal_node_handle; + typedef typename Kd_tree::Internal_node_const_handle Internal_node_const_handle; + typedef typename Kd_tree::Leaf_node_handle Leaf_node_handle; + typedef typename Kd_tree::Leaf_node_const_handle Leaf_node_const_handle; + typedef typename TreeTraits::Point_d Point_d; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree::Separator Separator; + typedef typename Kd_tree::Point_d_iterator Point_d_iterator; + typedef typename Kd_tree::iterator iterator; + typedef typename Kd_tree::D D; + + bool leaf; + + public : + Kd_tree_node(bool leaf_) + :leaf(leaf_){} + + bool is_leaf() const{ + return leaf; + } + + std::size_t + num_items() const + { + if (is_leaf()){ + Leaf_node_const_handle node = + static_cast(this); + return node->size(); + } + else { + Internal_node_const_handle node = + static_cast(this); + return node->lower()->num_items() + node->upper()->num_items(); + } + } + + std::size_t + num_nodes() const + { + if (is_leaf()) return 1; + else { + Internal_node_const_handle node = + static_cast(this); + return node->lower()->num_nodes() + node->upper()->num_nodes(); + } + } + + int + depth(const int current_max_depth) const + { + if (is_leaf()){ + return current_max_depth; + } + else { + Internal_node_const_handle node = + static_cast(this); + return + (std::max)( node->lower()->depth(current_max_depth + 1), + node->upper()->depth(current_max_depth + 1)); + } + } + + int + depth() const + { + return depth(1); + } + + template + OutputIterator + tree_items(OutputIterator it) const { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + {*it=*i; ++it;} + } + else { + Internal_node_const_handle node = + static_cast(this); + it=node->lower()->tree_items(it); + it=node->upper()->tree_items(it); + } + return it; + } + + + boost::optional + any_tree_item() const { + boost::optional result = boost::none; + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast(this); + if (node->size()>0){ + return boost::make_optional(*(node->begin())); + } + } + else { + Internal_node_const_handle node = + static_cast(this); + result = node->lower()->any_tree_item(); + if(! result){ + result = node->upper()->any_tree_item(); + } + } + return result; + } + + + void + indent(int d) const + { + for(int i = 0; i < d; i++){ + std::cout << " "; + } + } + + + void + print(int d = 0) const + { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast(this); + indent(d); + std::cout << "leaf" << std::endl; + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + {indent(d);std::cout << *i << std::endl;} + } + else { + Internal_node_const_handle node = + static_cast(this); + indent(d); + std::cout << "lower tree" << std::endl; + node->lower()->print(d+1); + indent(d); + std::cout << "separator: dim = " << node->cutting_dimension() << " val = " << node->cutting_value() << std::endl; + indent(d); + std::cout << "upper tree" << std::endl; + node->upper()->print(d+1); + } + } + + + template + OutputIterator + search(OutputIterator it, const FuzzyQueryItem& q, + Kd_tree_rectangle& b) const + { + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + if (q.contains(*i)) + {*it++=*i;} + } + else { + Internal_node_const_handle node = + static_cast(this); + // after splitting b denotes the lower part of b + Kd_tree_rectangle b_upper(b); + b.split(b_upper, node->cutting_dimension(), + node->cutting_value()); + + if (q.outer_range_contains(b)) + it=node->lower()->tree_items(it); + else + if (q.inner_range_intersects(b)) + it=node->lower()->search(it,q,b); + if (q.outer_range_contains(b_upper)) + it=node->upper()->tree_items(it); + else + if (q.inner_range_intersects(b_upper)) + it=node->upper()->search(it,q,b_upper); + }; + return it; + } + + + template + boost::optional + search_any_point(const FuzzyQueryItem& q, + Kd_tree_rectangle& b) const + { + boost::optional result = boost::none; + if (is_leaf()) { + Leaf_node_const_handle node = + static_cast(this); + if (node->size()>0) + for (iterator i=node->begin(); i != node->end(); i++) + if (q.contains(*i)) + { result = *i; break; } + } + else { + Internal_node_const_handle node = + static_cast(this); + // after splitting b denotes the lower part of b + Kd_tree_rectangle b_upper(b); + b.split(b_upper, node->cutting_dimension(), + node->cutting_value()); + + if (q.outer_range_contains(b)){ + result = node->lower()->any_tree_item(); + }else{ + if (q.inner_range_intersects(b)){ + result = node->lower()->search_any_point(q,b); + } + } + if(result){ + return result; + } + if (q.outer_range_contains(b_upper)){ + result = node->upper()->any_tree_item(); + }else{ + if (q.inner_range_intersects(b_upper)) + result = node->upper()->search_any_point(q,b_upper); + } + } + return result; + } + + }; + + + template < class TreeTraits, class Splitter, class UseExtendedNode > + class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ + + friend class Kd_tree; + + typedef typename Kd_tree::iterator iterator; + typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; + typedef typename TreeTraits::Point_d Point_d; + + private: + + // private variables for leaf nodes + boost::int32_t n; // denotes number of items in a leaf node + iterator data; // iterator to data in leaf node + + + public: + + // default constructor + Kd_tree_leaf_node() + {} + + Kd_tree_leaf_node(bool leaf_ ) + : Base(leaf_) + {} + + Kd_tree_leaf_node(bool leaf_,unsigned int n_ ) + : Base(leaf_), n(n_) + {} + + // members for all nodes + + // members for leaf nodes only + inline + unsigned int + size() const + { + return n; + } + + inline + iterator + begin() const + { + return data; + } + + inline + iterator + end() const + { + return data + n; + } + + inline + void + drop_last_point() + { + --n; + } + + }; //leaf node + + + + template < class TreeTraits, class Splitter, class UseExtendedNode> + class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ + + friend class Kd_tree; + + typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; + typedef typename Kd_tree::Node_handle Node_handle; + typedef typename Kd_tree::Node_const_handle Node_const_handle; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree::Separator Separator; + + private: + + // private variables for internal nodes + boost::int32_t cut_dim; + FT cut_val; + Node_handle lower_ch, upper_ch; + + + // private variables for extended internal nodes + FT upper_low_val; + FT upper_high_val; + FT lower_low_val; + FT lower_high_val; + + + public: + + // default constructor + Kd_tree_internal_node() + {} + + Kd_tree_internal_node(bool leaf_) + : Base(leaf_) + {} + + + // members for internal node and extended internal node + + inline + Node_const_handle + lower() const + { + return lower_ch; + } + + inline + Node_const_handle + upper() const + { + return upper_ch; + } + + inline + Node_handle + lower() + { + return lower_ch; + } + + inline + Node_handle + upper() + { + return upper_ch; + } + + inline + void + set_lower(Node_handle nh) + { + lower_ch = nh; + } + + inline + void + set_upper(Node_handle nh) + { + upper_ch = nh; + } + + // inline Separator& separator() {return sep; } + // use instead + inline + void set_separator(Separator& sep){ + cut_dim = sep.cutting_dimension(); + cut_val = sep.cutting_value(); + } + + inline + FT + cutting_value() const + { + return cut_val; + } + + inline + int + cutting_dimension() const + { + return cut_dim; + } + + // members for extended internal node only + inline + FT + upper_low_value() const + { + return upper_low_val; + } + + inline + FT + upper_high_value() const + { + return upper_high_val; + } + + inline + FT + lower_low_value() const + { + return lower_low_val; + } + + inline + FT + lower_high_value() const + { + return lower_high_val; + } + + /*Separator& + separator() + { + return Separator(cutting_dimension,cutting_value); + }*/ + + + };//internal node + + template < class TreeTraits, class Splitter> + class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, Tag_false >{ + + friend class Kd_tree; + + typedef Kd_tree_node< TreeTraits, Splitter, Tag_false> Base; + typedef typename Kd_tree::Node_handle Node_handle; + typedef typename Kd_tree::Node_const_handle Node_const_handle; + + typedef typename TreeTraits::FT FT; + typedef typename Kd_tree::Separator Separator; + + private: + + // private variables for internal nodes + boost::uint8_t cut_dim; + FT cut_val; + + Node_handle lower_ch, upper_ch; + + public: + + // default constructor + Kd_tree_internal_node() + {} + + Kd_tree_internal_node(bool leaf_) + : Base(leaf_) + {} + + + // members for internal node and extended internal node + + inline + Node_const_handle + lower() const + { + return lower_ch; + } + + inline + Node_const_handle + upper() const + { + return upper_ch; + } + + inline + Node_handle + lower() + { + return lower_ch; + } + + inline + Node_handle + upper() + { + return upper_ch; + } + + inline + void + set_lower(Node_handle nh) + { + lower_ch = nh; + } + + inline + void + set_upper(Node_handle nh) + { + upper_ch = nh; + } + + // inline Separator& separator() {return sep; } + // use instead + + inline + void set_separator(Separator& sep){ + cut_dim = sep.cutting_dimension(); + cut_val = sep.cutting_value(); + } + + inline + FT + cutting_value() const + { + return cut_val; + } + + inline + int + cutting_dimension() const + { + return cut_dim; + } + + /* Separator& + separator() + { + return Separator(cutting_dimension,cutting_value); + }*/ + + + };//internal node + + + +} // namespace CGAL +#endif // CGAL_KDTREE_NODE_H diff --git a/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h b/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h new file mode 100644 index 00000000..dbe707ed --- /dev/null +++ b/src/common/include/gudhi_patches/CGAL/Orthogonal_incremental_neighbor_search.h @@ -0,0 +1,621 @@ +// Copyright (c) 2002,2011 Utrecht University (The Netherlands). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Hans Tangelder () + +#ifndef CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH +#define CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH + +#include "CGAL/Kd_tree.h" + +#include +#include +#include +#include +#include +#include + +namespace CGAL { + + template ::type, + class Splitter_ = Sliding_midpoint, + class Tree_= Kd_tree > + class Orthogonal_incremental_neighbor_search { + + public: + typedef Splitter_ Splitter; + typedef Tree_ Tree; + typedef Distance_ Distance; + typedef typename SearchTraits::Point_d Point_d; + typedef typename Distance::Query_item Query_item; + typedef typename SearchTraits::FT FT; + typedef typename Tree::Point_d_iterator Point_d_iterator; + typedef typename Tree::Node_const_handle Node_const_handle; + + typedef std::pair Point_with_transformed_distance; + typedef CGAL::cpp11::tuple > Node_with_distance; + typedef std::vector Node_with_distance_vector; + typedef std::vector Point_with_transformed_distance_vector; + + template + struct Object_wrapper + { + T object; + Object_wrapper(const T& t):object(t){} + const T& operator* () const { return object; } + const T* operator-> () const { return &object; } + }; + + class Iterator_implementation { + SearchTraits traits; + public: + + int number_of_neighbours_computed; + int number_of_internal_nodes_visited; + int number_of_leaf_nodes_visited; + int number_of_items_visited; + + private: + + typedef std::vector Distance_vector; + + Distance_vector dists; + + Distance Orthogonal_distance_instance; + + FT multiplication_factor; + + Query_item query_point; + + FT distance_to_root; + + bool search_nearest_neighbour; + + FT rd; + + + class Priority_higher { + public: + + bool search_nearest; + + Priority_higher(bool search_the_nearest_neighbour) + : search_nearest(search_the_nearest_neighbour) + {} + + //highest priority is smallest distance + bool + operator() (Node_with_distance* n1, Node_with_distance* n2) const + { + return (search_nearest) ? (CGAL::cpp11::get<1>(*n1) > CGAL::cpp11::get<1>(*n2)) : (CGAL::cpp11::get<1>(*n2) > CGAL::cpp11::get<1>(*n1)); + } + }; + + class Distance_smaller { + + public: + + bool search_nearest; + + Distance_smaller(bool search_the_nearest_neighbour) + : search_nearest(search_the_nearest_neighbour) + {} + + //highest priority is smallest distance + bool operator() (Point_with_transformed_distance* p1, Point_with_transformed_distance* p2) const + { + return (search_nearest) ? (p1->second > p2->second) : (p2->second > p1->second); + } + }; + + + std::priority_queue PriorityQueue; + + public: + std::priority_queue Item_PriorityQueue; + + + public: + + int reference_count; + + + + // constructor + Iterator_implementation(const Tree& tree,const Query_item& q, const Distance& tr, + FT Eps=FT(0.0), bool search_nearest=true) + : traits(tree.traits()),number_of_neighbours_computed(0), number_of_internal_nodes_visited(0), + number_of_leaf_nodes_visited(0), number_of_items_visited(0), + Orthogonal_distance_instance(tr), multiplication_factor(Orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)), + query_point(q), search_nearest_neighbour(search_nearest), + PriorityQueue(Priority_higher(search_nearest)), Item_PriorityQueue(Distance_smaller(search_nearest)), + reference_count(1) + + + { + if (tree.empty()) return; + + typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits.construct_cartesian_const_iterator_d_object(); + int dim = static_cast(std::distance(ccci(q), ccci(q,0))); + + dists.resize(dim); + for(int i=0 ; i(*The_Root); + Compute_the_next_nearest_neighbour(); + } + else{ + distance_to_root= + Orthogonal_distance_instance.max_distance_to_rectangle(q, + tree.bounding_box(), dists); + Node_with_distance *The_Root = new Node_with_distance(tree.root(), + distance_to_root, dists); + PriorityQueue.push(The_Root); + + // rd is the distance of the top of the priority queue to q + rd=CGAL::cpp11::get<1>(*The_Root); + Compute_the_next_furthest_neighbour(); + } + + + } + + // * operator + const Point_with_transformed_distance& + operator* () const + { + return *(Item_PriorityQueue.top()); + } + + // prefix operator + Iterator_implementation& + operator++() + { + Delete_the_current_item_top(); + if(search_nearest_neighbour) + Compute_the_next_nearest_neighbour(); + else + Compute_the_next_furthest_neighbour(); + return *this; + } + + // postfix operator + Object_wrapper + operator++(int) + { + Object_wrapper result( *(Item_PriorityQueue.top()) ); + ++*this; + return result; + } + + // Print statistics of the general priority search process. + std::ostream& + statistics (std::ostream& s) const { + s << "Orthogonal priority search statistics:" + << std::endl; + s << "Number of internal nodes visited:" + << number_of_internal_nodes_visited << std::endl; + s << "Number of leaf nodes visited:" + << number_of_leaf_nodes_visited << std::endl; + s << "Number of items visited:" + << number_of_items_visited << std::endl; + s << "Number of neighbours computed:" + << number_of_neighbours_computed << std::endl; + return s; + } + + + //destructor + ~Iterator_implementation() + { + while (!PriorityQueue.empty()) { + Node_with_distance* The_top=PriorityQueue.top(); + PriorityQueue.pop(); + delete The_top; + } + while (!Item_PriorityQueue.empty()) { + Point_with_transformed_distance* The_top=Item_PriorityQueue.top(); + Item_PriorityQueue.pop(); + delete The_top; + } + } + + private: + + void + Delete_the_current_item_top() + { + Point_with_transformed_distance* The_item_top=Item_PriorityQueue.top(); + Item_PriorityQueue.pop(); + delete The_item_top; + } + + void + Compute_the_next_nearest_neighbour() + { + // compute the next item + bool next_neighbour_found=false; + if (!(Item_PriorityQueue.empty())) { + next_neighbour_found= + (multiplication_factor*rd > Item_PriorityQueue.top()->second); + } + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); + // otherwise browse the tree further + while ((!next_neighbour_found) && (!PriorityQueue.empty())) { + Node_with_distance* The_node_top=PriorityQueue.top(); + Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); + dists = CGAL::cpp11::get<2>(*The_node_top); + PriorityQueue.pop(); + delete The_node_top; + FT copy_rd=rd; + while (!(N->is_leaf())) { // compute new distance + typename Tree::Internal_node_const_handle node = + static_cast(N); + number_of_internal_nodes_visited++; + int new_cut_dim=node->cutting_dimension(); + FT new_rd,dst = dists[new_cut_dim]; + FT val = *(query_point_it + new_cut_dim); + FT diff1 = val - node->upper_low_value(); + FT diff2 = val - node->lower_high_value(); + if (diff1 + diff2 < FT(0.0)) { + new_rd= + Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); + + CGAL_assertion(new_rd >= copy_rd); + dists[new_cut_dim] = diff1; + Node_with_distance *Upper_Child = + new Node_with_distance(node->upper(), new_rd, dists); + PriorityQueue.push(Upper_Child); + dists[new_cut_dim] = dst; + N=node->lower(); + + } + else { // compute new distance + new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); + CGAL_assertion(new_rd >= copy_rd); + dists[new_cut_dim] = diff2; + Node_with_distance *Lower_Child = + new Node_with_distance(node->lower(), new_rd, dists); + PriorityQueue.push(Lower_Child); + dists[new_cut_dim] = dst; + N=node->upper(); + } + } + // n is a leaf + typename Tree::Leaf_node_const_handle node = + static_cast(N); + number_of_leaf_nodes_visited++; + if (node->size() > 0) { + for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { + number_of_items_visited++; + FT distance_to_query_point= + Orthogonal_distance_instance.transformed_distance(query_point,*it); + Point_with_transformed_distance *NN_Candidate= + new Point_with_transformed_distance(*it,distance_to_query_point); + Item_PriorityQueue.push(NN_Candidate); + } + // old top of PriorityQueue has been processed, + // hence update rd + + if (!(PriorityQueue.empty())) { + rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); + next_neighbour_found = + (multiplication_factor*rd > + Item_PriorityQueue.top()->second); + } + else // priority queue empty => last neighbour found + { + next_neighbour_found=true; + } + + number_of_neighbours_computed++; + } + } // next_neighbour_found or priority queue is empty + // in the latter case also the item priority quee is empty + } + + + void + Compute_the_next_furthest_neighbour() + { + // compute the next item + bool next_neighbour_found=false; + if (!(Item_PriorityQueue.empty())) { + next_neighbour_found= + (rd < multiplication_factor*Item_PriorityQueue.top()->second); + } + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d query_point_it = construct_it(query_point); + // otherwise browse the tree further + while ((!next_neighbour_found) && (!PriorityQueue.empty())) { + Node_with_distance* The_node_top=PriorityQueue.top(); + Node_const_handle N= CGAL::cpp11::get<0>(*The_node_top); + dists = CGAL::cpp11::get<2>(*The_node_top); + PriorityQueue.pop(); + delete The_node_top; + FT copy_rd=rd; + while (!(N->is_leaf())) { // compute new distance + typename Tree::Internal_node_const_handle node = + static_cast(N); + number_of_internal_nodes_visited++; + int new_cut_dim=node->cutting_dimension(); + FT new_rd,dst = dists[new_cut_dim]; + FT val = *(query_point_it + new_cut_dim); + FT diff1 = val - node->upper_low_value(); + FT diff2 = val - node->lower_high_value(); + if (diff1 + diff2 < FT(0.0)) { + diff1 = val - node->upper_high_value(); + new_rd= + Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); + Node_with_distance *Lower_Child = + new Node_with_distance(node->lower(), copy_rd, dists); + PriorityQueue.push(Lower_Child); + N=node->upper(); + dists[new_cut_dim] = diff1; + copy_rd=new_rd; + + } + else { // compute new distance + diff2 = val - node->lower_low_value(); + new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); + Node_with_distance *Upper_Child = + new Node_with_distance(node->upper(), copy_rd, dists); + PriorityQueue.push(Upper_Child); + N=node->lower(); + dists[new_cut_dim] = diff2; + copy_rd=new_rd; + } + } + // n is a leaf + typename Tree::Leaf_node_const_handle node = + static_cast(N); + number_of_leaf_nodes_visited++; + if (node->size() > 0) { + for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { + number_of_items_visited++; + FT distance_to_query_point= + Orthogonal_distance_instance.transformed_distance(query_point,*it); + Point_with_transformed_distance *NN_Candidate= + new Point_with_transformed_distance(*it,distance_to_query_point); + Item_PriorityQueue.push(NN_Candidate); + } + // old top of PriorityQueue has been processed, + // hence update rd + + if (!(PriorityQueue.empty())) { + rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); + next_neighbour_found = + (multiplication_factor*rd < + Item_PriorityQueue.top()->second); + } + else // priority queue empty => last neighbour found + { + next_neighbour_found=true; + } + + number_of_neighbours_computed++; + } + } // next_neighbour_found or priority queue is empty + // in the latter case also the item priority quee is empty + } + }; // class Iterator_implementaion + + + + + + + + + + public: + class iterator; + typedef iterator const_iterator; + + // constructor + Orthogonal_incremental_neighbor_search(const Tree& tree, + const Query_item& q, FT Eps = FT(0.0), + bool search_nearest=true, const Distance& tr=Distance()) + : m_tree(tree),m_query(q),m_dist(tr),m_Eps(Eps),m_search_nearest(search_nearest) + {} + + iterator + begin() const + { + return iterator(m_tree,m_query,m_dist,m_Eps,m_search_nearest); + } + + iterator + end() const + { + return iterator(); + } + + std::ostream& + statistics(std::ostream& s) + { + begin()->statistics(s); + return s; + } + + + + + class iterator { + + public: + + typedef std::input_iterator_tag iterator_category; + typedef Point_with_transformed_distance value_type; + typedef Point_with_transformed_distance* pointer; + typedef const Point_with_transformed_distance& reference; + typedef std::size_t size_type; + typedef std::ptrdiff_t difference_type; + typedef int distance_type; + + //class Iterator_implementation; + Iterator_implementation *Ptr_implementation; + + + public: + + // default constructor + iterator() + : Ptr_implementation(0) + {} + + int + the_number_of_items_visited() + { + return Ptr_implementation->number_of_items_visited; + } + + // constructor + iterator(const Tree& tree,const Query_item& q, const Distance& tr=Distance(), FT eps=FT(0.0), + bool search_nearest=true) + : Ptr_implementation(new Iterator_implementation(tree, q, tr, eps, search_nearest)) + {} + + // copy constructor + iterator(const iterator& Iter) + { + Ptr_implementation = Iter.Ptr_implementation; + if (Ptr_implementation != 0) Ptr_implementation->reference_count++; + } + + iterator& operator=(const iterator& Iter) + { + if (Ptr_implementation != Iter.Ptr_implementation){ + if (Ptr_implementation != 0 && --(Ptr_implementation->reference_count)==0) { + delete Ptr_implementation; + } + Ptr_implementation = Iter.Ptr_implementation; + if (Ptr_implementation != 0) Ptr_implementation->reference_count++; + } + return *this; + } + + + const Point_with_transformed_distance& + operator* () const + { + return *(*Ptr_implementation); + } + + // -> operator + const Point_with_transformed_distance* + operator-> () const + { + return &*(*Ptr_implementation); + } + + // prefix operator + iterator& + operator++() + { + ++(*Ptr_implementation); + return *this; + } + + // postfix operator + Object_wrapper + operator++(int) + { + return (*Ptr_implementation)++; + } + + + bool + operator==(const iterator& It) const + { + if ( + ((Ptr_implementation == 0) || + Ptr_implementation->Item_PriorityQueue.empty()) && + ((It.Ptr_implementation == 0) || + It.Ptr_implementation->Item_PriorityQueue.empty()) + ) + return true; + // else + return (Ptr_implementation == It.Ptr_implementation); + } + + bool + operator!=(const iterator& It) const + { + return !(*this == It); + } + + std::ostream& + statistics (std::ostream& s) + { + Ptr_implementation->statistics(s); + return s; + } + + ~iterator() + { + if (Ptr_implementation != 0) { + Ptr_implementation->reference_count--; + if (Ptr_implementation->reference_count==0) { + delete Ptr_implementation; + Ptr_implementation = 0; + } + } + } + + + }; // class iterator + + //data members + const Tree& m_tree; + Query_item m_query; + Distance m_dist; + FT m_Eps; + bool m_search_nearest; + }; // class + + template + void swap (typename Orthogonal_incremental_neighbor_search::iterator& x, + typename Orthogonal_incremental_neighbor_search::iterator& y) + { + typename Orthogonal_incremental_neighbor_search::iterator::Iterator_implementation + *tmp = x.Ptr_implementation; + x.Ptr_implementation = y.Ptr_implementation; + y.Ptr_implementation = tmp; + } + +} // namespace CGAL + +#endif // CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH_H diff --git a/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt b/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt new file mode 100644 index 00000000..5b9581a0 --- /dev/null +++ b/src/common/include/gudhi_patches/Tangential_complex_CGAL_patches.txt @@ -0,0 +1,82 @@ +CGAL/Regular_triangulation_traits_adapter.h +CGAL/Triangulation_ds_vertex.h +CGAL/Triangulation_data_structure.h +CGAL/transforming_pair_iterator.h +CGAL/NewKernel_d/static_int.h +CGAL/NewKernel_d/Cartesian_LA_functors.h +CGAL/NewKernel_d/Cartesian_change_FT.h +CGAL/NewKernel_d/Wrapper/Vector_d.h +CGAL/NewKernel_d/Wrapper/Hyperplane_d.h +CGAL/NewKernel_d/Wrapper/Ref_count_obj.h +CGAL/NewKernel_d/Wrapper/Cartesian_wrap.h +CGAL/NewKernel_d/Wrapper/Point_d.h +CGAL/NewKernel_d/Wrapper/Segment_d.h +CGAL/NewKernel_d/Wrapper/Weighted_point_d.h +CGAL/NewKernel_d/Wrapper/Sphere_d.h +CGAL/NewKernel_d/Cartesian_per_dimension.h +CGAL/NewKernel_d/Kernel_object_converter.h +CGAL/NewKernel_d/KernelD_converter.h +CGAL/NewKernel_d/Vector/sse2.h +CGAL/NewKernel_d/Vector/avx4.h +CGAL/NewKernel_d/Vector/determinant_of_vectors_small_dim_internal.h +CGAL/NewKernel_d/Vector/determinant_of_iterator_to_points_from_points.h +CGAL/NewKernel_d/Vector/determinant_of_points_from_vectors.h +CGAL/NewKernel_d/Vector/array.h +CGAL/NewKernel_d/Vector/determinant_of_iterator_to_points_from_iterator_to_vectors.h +CGAL/NewKernel_d/Vector/determinant_of_iterator_to_vectors_from_vectors.h +CGAL/NewKernel_d/Vector/determinant_of_vectors_small_dim.h +CGAL/NewKernel_d/Vector/vector.h +CGAL/NewKernel_d/Vector/v2int.h +CGAL/NewKernel_d/Vector/mix.h +CGAL/NewKernel_d/Cartesian_static_filters.h +CGAL/NewKernel_d/Cartesian_LA_base.h +CGAL/NewKernel_d/Lazy_cartesian.h +CGAL/NewKernel_d/Coaffine.h +CGAL/NewKernel_d/store_kernel.h +CGAL/NewKernel_d/Dimension_base.h +CGAL/NewKernel_d/Kernel_3_interface.h +CGAL/NewKernel_d/Cartesian_complete.h +CGAL/NewKernel_d/Cartesian_base.h +CGAL/NewKernel_d/Cartesian_filter_K.h +CGAL/NewKernel_d/functor_tags.h +CGAL/NewKernel_d/Filtered_predicate2.h +CGAL/NewKernel_d/functor_properties.h +CGAL/NewKernel_d/Define_kernel_types.h +CGAL/NewKernel_d/LA_eigen/LA.h +CGAL/NewKernel_d/LA_eigen/constructors.h +CGAL/NewKernel_d/Types/Aff_transformation.h +CGAL/NewKernel_d/Types/Sphere.h +CGAL/NewKernel_d/Types/Hyperplane.h +CGAL/NewKernel_d/Types/Line.h +CGAL/NewKernel_d/Types/Ray.h +CGAL/NewKernel_d/Types/Iso_box.h +CGAL/NewKernel_d/Types/Weighted_point.h +CGAL/NewKernel_d/Types/Segment.h +CGAL/NewKernel_d/Kernel_d_interface.h +CGAL/NewKernel_d/utils.h +CGAL/NewKernel_d/Kernel_2_interface.h +CGAL/NewKernel_d/Cartesian_filter_NT.h +CGAL/NewKernel_d/function_objects_cartesian.h +CGAL/Convex_hull.h +CGAL/Triangulation_ds_full_cell.h +CGAL/Regular_triangulation.h +CGAL/Epick_d.h +CGAL/transforming_iterator.h +CGAL/iterator_from_indices.h +CGAL/Delaunay_triangulation.h +CGAL/IO/Triangulation_off_ostream.h +CGAL/typeset.h +CGAL/Triangulation_full_cell.h +CGAL/Triangulation.h +CGAL/internal/Static_or_dynamic_array.h +CGAL/internal/Combination_enumerator.h +CGAL/internal/Triangulation/utilities.h +CGAL/internal/Triangulation/Triangulation_ds_iterators.h +CGAL/internal/Triangulation/Dummy_TDS.h +CGAL/argument_swaps.h +CGAL/Epeck_d.h +CGAL/determinant_of_vectors.h +CGAL/TDS_full_cell_default_storage_policy.h +CGAL/TDS_full_cell_mirror_storage_policy.h +CGAL/Triangulation_face.h +CGAL/Triangulation_vertex.h -- cgit v1.2.3 From a823bfcb70ed76e8858604050570ff8fe33f6667 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Fri, 20 Jan 2017 12:42:27 +0000 Subject: cpplint fixes doxygen warning fixes git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/bottleneck_misc_fixes@1972 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 5b2328016edb6b5b30de668a0488c576d3d92c40 --- data/points/generator/hypergenerator.cpp | 16 +++++++--------- .../example/bottleneck_basic_example.cpp | 6 ++++-- .../example/bottleneck_read_file_example.cpp | 2 ++ src/Bottleneck_distance/include/gudhi/Bottleneck.h | 9 +++++++-- src/Bottleneck_distance/include/gudhi/Graph_matching.h | 3 +++ .../include/gudhi/Neighbors_finder.h | 1 + .../include/gudhi/Persistence_graph.h | 10 +++++++--- src/Doxyfile | 3 ++- src/common/doc/main_page.h | 17 +++++++++++------ .../example/example_CGAL_3D_points_off_reader.cpp | 2 +- src/common/example/example_CGAL_points_off_reader.cpp | 2 +- 11 files changed, 46 insertions(+), 25 deletions(-) (limited to 'src/Bottleneck_distance/include/gudhi/Bottleneck.h') diff --git a/data/points/generator/hypergenerator.cpp b/data/points/generator/hypergenerator.cpp index 60890b44..5831de18 100644 --- a/data/points/generator/hypergenerator.cpp +++ b/data/points/generator/hypergenerator.cpp @@ -29,7 +29,7 @@ #include #include #include // for std::ofstream - +#include typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > K; typedef K::Point_d Point; @@ -47,24 +47,22 @@ int main(int argc, char **argv) { usage(argv[0]); } - int points_number = 0; - int returnedScanValue = sscanf(argv[4], "%d", &points_number); - if ((returnedScanValue == EOF) || (points_number <= 0)) { + int points_number = atoi(argv[4]); + if (points_number <= 0) { std::cerr << "Error: " << argv[4] << " is not correct" << std::endl; usage(argv[0]); } - int dimension = 0; - returnedScanValue = sscanf(argv[5], "%d", &dimension); - if ((returnedScanValue == EOF) || (dimension <= 0)) { + int dimension = atoi(argv[5]); + if (dimension <= 0) { std::cerr << "Error: " << argv[5] << " is not correct" << std::endl; usage(argv[0]); } double radius = 1.0; if (argc == 7) { - returnedScanValue = sscanf(argv[6], "%lf", &radius); - if ((returnedScanValue == EOF) || (radius <= 0.0)) { + radius = atof(argv[6]); + if (radius <= 0.0) { std::cerr << "Error: " << argv[6] << " is not correct" << std::endl; usage(argv[0]); } diff --git a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp index 91a7302f..d0ca4e20 100644 --- a/src/Bottleneck_distance/example/bottleneck_basic_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_basic_example.cpp @@ -21,10 +21,13 @@ */ #include + #include +#include +#include // for pair +#include // for numeric_limits int main() { - std::vector< std::pair > v1, v2; v1.emplace_back(2.7, 3.7); @@ -44,5 +47,4 @@ int main() { b = Gudhi::persistence_diagram::bottleneck_distance(v1, v2, 0.1); std::cout << "Approx bottleneck distance = " << b << std::endl; - } diff --git a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp index 4c74b66e..bde05825 100644 --- a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp +++ b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp @@ -24,6 +24,8 @@ #include #include +#include +#include // for pair #include #include #include diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h index 2b7e4767..b5641e29 100644 --- a/src/Bottleneck_distance/include/gudhi/Bottleneck.h +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -24,6 +24,11 @@ #define BOTTLENECK_H_ #include + +#include +#include // for max +#include // for numeric_limits + #include namespace Gudhi { @@ -41,7 +46,7 @@ double bottleneck_distance_approx(Persistence_graph& g, double e) { if (step <= b_lower_bound || step >= b_upper_bound) // Avoid precision problem break; m.set_r(step); - while (m.multi_augment()); // compute a maximum matching (in the graph corresponding to the current r) + while (m.multi_augment()) {}; // compute a maximum matching (in the graph corresponding to the current r) if (m.perfect()) { m = biggest_unperfect; b_upper_bound = step; @@ -63,7 +68,7 @@ double bottleneck_distance_exact(Persistence_graph& 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) + 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; diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h index 253c89b4..e1708c5b 100644 --- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -25,6 +25,9 @@ #include +#include +#include + namespace Gudhi { namespace persistence_diagram { diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h index 96ece360..cd5486f8 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -34,6 +34,7 @@ #include #include +#include namespace Gudhi { diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h index 3a4a5fec..39efc082 100644 --- a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h +++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h @@ -23,9 +23,11 @@ #ifndef PERSISTENCE_GRAPH_H_ #define PERSISTENCE_GRAPH_H_ +#include + #include #include -#include +#include // for numeric_limits namespace Gudhi { @@ -92,10 +94,12 @@ Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1, swap(u, v); std::sort(u_alive.begin(), u_alive.end()); std::sort(v_alive.begin(), v_alive.end()); - if (u_alive.size() != v_alive.size()) + 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)); + } } inline bool Persistence_graph::on_the_u_diagonal(int u_point_index) const { diff --git a/src/Doxyfile b/src/Doxyfile index 71aa038d..bf7deb57 100644 --- a/src/Doxyfile +++ b/src/Doxyfile @@ -784,7 +784,8 @@ EXCLUDE = data/ \ example/ \ GudhUI/ \ cmake/ \ - debian/ + debian/ \ + include/gudhi_patches/ # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or # directories that are symbolic links (a Unix file system feature) are excluded diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h index cf67b558..660327d9 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -261,9 +261,11 @@ make \endverbatim * Persistent_cohomology/custom_persistence_sort.cpp * * The following example requires CGAL version ≥ 4.8.0: - * \li - * - * Bottleneck_distance/bottleneck_example.cpp + * \li + * Bottleneck_distance/bottleneck_basic_example.cpp + * \li + * Bottleneck_distance/bottleneck_read_file_example.cpp + * * \subsection eigen3 Eigen3 * The \ref alpha_complex data structure and few examples requires * Eigen3 is a C++ template library for linear algebra: @@ -275,8 +277,10 @@ make \endverbatim * Alpha_complex/Alpha_complex_from_off.cpp * \li * Alpha_complex/Alpha_complex_from_points.cpp - * \li - * Bottleneck_distance/bottleneck_example.cpp + * \li + * Bottleneck_distance/bottleneck_basic_example.cpp + * \li + * Bottleneck_distance/bottleneck_read_file_example.cpp * \li * Persistent_cohomology/alpha_complex_persistence.cpp * \li @@ -358,7 +362,8 @@ make \endverbatim /*! @file Examples * @example Alpha_complex/Alpha_complex_from_off.cpp * @example Alpha_complex/Alpha_complex_from_points.cpp - * @example Bottleneck_distance/bottleneck_example.cpp + * @example Bottleneck_distance/bottleneck_basic_example.cpp + * @example Bottleneck_distance/bottleneck_read_file_example.cpp * @example Bitmap_cubical_complex/Bitmap_cubical_complex.cpp * @example Bitmap_cubical_complex/Bitmap_cubical_complex_periodic_boundary_conditions.cpp * @example Bitmap_cubical_complex/Random_bitmap_cubical_complex.cpp diff --git a/src/common/example/example_CGAL_3D_points_off_reader.cpp b/src/common/example/example_CGAL_3D_points_off_reader.cpp index d48bb17d..665b7a29 100644 --- a/src/common/example/example_CGAL_3D_points_off_reader.cpp +++ b/src/common/example/example_CGAL_3D_points_off_reader.cpp @@ -32,7 +32,7 @@ int main(int argc, char **argv) { // Retrieve the triangulation std::vector point_cloud = off_reader.get_point_cloud(); - int n {0}; + int n {}; for (auto point : point_cloud) { ++n; std::cout << "Point[" << n << "] = (" << point[0] << ", " << point[1] << ", " << point[2] << ")\n"; diff --git a/src/common/example/example_CGAL_points_off_reader.cpp b/src/common/example/example_CGAL_points_off_reader.cpp index 4522174a..8c6a6b54 100644 --- a/src/common/example/example_CGAL_points_off_reader.cpp +++ b/src/common/example/example_CGAL_points_off_reader.cpp @@ -34,7 +34,7 @@ int main(int argc, char **argv) { // Retrieve the triangulation std::vector point_cloud = off_reader.get_point_cloud(); - int n {0}; + int n {}; for (auto point : point_cloud) { std::cout << "Point[" << n << "] = "; for (std::size_t i {0}; i < point.size(); i++) -- cgit v1.2.3