diff options
author | vrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2017-04-14 05:49:23 +0000 |
---|---|---|
committer | vrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2017-04-14 05:49:23 +0000 |
commit | 30bbfd7d261b098efdca0fe16f4968040effd7ed (patch) | |
tree | 667bbb1e4c37fedbac5547296fca4007a9a825cf /src | |
parent | d9922fab02fd28d60601dab1aeb5ede2246286d7 (diff) | |
parent | 387a5af3ee1b664346eb9686f00c986e9f7a1e3e (diff) |
Merge last trunk modifications
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/ST_cythonize@2346 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: 29b4678c6e1ecaec8c43fc68a36aa988ddcbf298
Diffstat (limited to 'src')
22 files changed, 198 insertions, 142 deletions
diff --git a/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp b/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp index fd9f0858..fd164b22 100644 --- a/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp +++ b/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp @@ -74,7 +74,7 @@ int main(int argc, char * argv[]) { // -------------------------------------------- // Rips persistence // -------------------------------------------- - Rips_complex rips_complex(off_reader.get_point_cloud(), threshold, Euclidean_distance()); + Rips_complex rips_complex(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); // Construct the Rips complex in a Simplex Tree Simplex_tree rips_stree; diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h index e1708c5b..f51e22e9 100644 --- a/src/Bottleneck_distance/include/gudhi/Graph_matching.h +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -26,7 +26,8 @@ #include <gudhi/Neighbors_finder.h> #include <vector> -#include <list> +#include <unordered_set> +#include <algorithm> namespace Gudhi { @@ -40,8 +41,6 @@ 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. */ @@ -50,12 +49,12 @@ class Graph_matching { void set_r(double r); private: - Persistence_graph& g; + Persistence_graph* gp; double r; /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */ std::vector<int> v_to_u; /** \internal \brief All the unmatched points in U. */ - std::list<int> unmatched_in_u; + std::unordered_set<int> unmatched_in_u; /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ Layered_neighbors_finder layering() const; @@ -66,17 +65,9 @@ class Graph_matching { }; inline Graph_matching::Graph_matching(Persistence_graph& g) - : g(g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u() { + : gp(&g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u(g.size()) { 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; + unmatched_in_u.insert(u_point_index); } inline bool Graph_matching::perfect() const { @@ -88,12 +79,12 @@ 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(gp->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<int> tries(unmatched_in_u); + std::vector<int> tries(unmatched_in_u.cbegin(), unmatched_in_u.cend()); 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; @@ -133,12 +124,12 @@ inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u } inline Layered_neighbors_finder Graph_matching::layering() const { - std::list<int> u_vertices(unmatched_in_u); - std::list<int> v_vertices; - Neighbors_finder nf(g, r); - for (int v_point_index = 0; v_point_index < g.size(); ++v_point_index) + std::vector<int> u_vertices(unmatched_in_u.cbegin(), unmatched_in_u.cend()); + std::vector<int> v_vertices; + Neighbors_finder nf(*gp, r); + for (int v_point_index = 0; v_point_index < gp->size(); ++v_point_index) nf.add(v_point_index); - Layered_neighbors_finder layered_nf(g, r); + Layered_neighbors_finder layered_nf(*gp, 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) { @@ -166,7 +157,8 @@ inline Layered_neighbors_finder Graph_matching::layering() const { } inline void Graph_matching::update(std::vector<int>& path) { - unmatched_in_u.remove(path.front()); + // Must return 1. + unmatched_in_u.erase(path.front()); for (auto it = path.cbegin(); it != path.cend(); ++it) { // Be careful, the iterator is incremented twice each time int tmp = *it; diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h index cd5486f8..bdc47578 100644 --- a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -25,9 +25,6 @@ // Inclusion order is important for CGAL patch #include <CGAL/Kd_tree.h> -#include <CGAL/Kd_tree_node.h> -#include <CGAL/Orthogonal_k_neighbor_search.h> -#include <CGAL/Weighted_Minkowski_distance.h> #include <CGAL/Search_traits.h> #include <gudhi/Persistence_graph.h> @@ -40,6 +37,33 @@ namespace Gudhi { namespace persistence_diagram { +/** \internal \brief Variant of CGAL::Fuzzy_iso_box to ensure that the box ic closed. It isn't clear how necessary that is. + */ +struct Square_query { + typedef CGAL::Dimension_tag<2> D; + typedef Internal_point Point_d; + typedef double FT; + bool contains(Point_d p) const { + return std::abs(p.x()-c.x())<=size && std::abs(p.y()-c.y())<=size; + } + bool inner_range_intersects(CGAL::Kd_tree_rectangle<FT,D> const&r) const { + return + r.max_coord(0) >= c.x() - size && + r.min_coord(0) <= c.x() + size && + r.max_coord(1) >= c.y() - size && + r.min_coord(1) <= c.y() + size; + } + bool outer_range_contains(CGAL::Kd_tree_rectangle<FT,D> const&r) const { + return + r.min_coord(0) >= c.x() - size && + r.max_coord(0) <= c.x() + size && + r.min_coord(1) >= c.y() - size && + r.max_coord(1) <= c.y() + size; + } + Point_d c; + FT size; +}; + /** \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). * @@ -51,9 +75,7 @@ namespace persistence_diagram { class Neighbors_finder { typedef CGAL::Dimension_tag<2> D; typedef CGAL::Search_traits<double, Internal_point, const double*, Construct_coord_iterator, D> Traits; - typedef CGAL::Weighted_Minkowski_distance<Traits> Distance; - typedef CGAL::Orthogonal_k_neighbor_search<Traits, Distance> K_neighbor_search; - typedef K_neighbor_search::Tree Kd_tree; + typedef CGAL::Kd_tree<Traits> Kd_tree; public: /** \internal \brief Constructor taking the near distance definition as parameter. */ @@ -123,15 +145,13 @@ 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::array<double, 2> 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) + auto neighbor = kd_t.search_any_point(Square_query{u_point, r}); + if(!neighbor) return null_point_index(); - tmp = it->first.point_index; - kd_t.remove(g.get_v_point(tmp)); + tmp = neighbor->point_index; + auto point = g.get_v_point(tmp); + int idx = point.point_index; + kd_t.remove(point, [idx](Internal_point const&p){return p.point_index == idx;}); } return tmp; } diff --git a/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp b/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp index ba752999..252e8aef 100644 --- a/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp +++ b/src/Persistent_cohomology/benchmark/performance_rips_persistence.cpp @@ -83,7 +83,7 @@ int main(int argc, char * argv[]) { // Compute the proximity graph of the points start = std::chrono::system_clock::now(); - Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Euclidean_distance()); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); end = std::chrono::system_clock::now(); elapsed_sec = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count(); std::cout << "Compute Rips graph in " << elapsed_sec << " ms.\n"; diff --git a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp index 7674b5a5..dae36ed2 100644 --- a/src/Persistent_cohomology/example/rips_multifield_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_multifield_persistence.cpp @@ -62,7 +62,7 @@ int main(int argc, char * argv[]) { program_options(argc, argv, off_file_points, filediag, threshold, dim_max, min_p, max_p, min_persistence); Points_off_reader off_reader(off_file_points); - Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Euclidean_distance()); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); // Construct the Rips complex in a Simplex Tree Simplex_tree simplex_tree; diff --git a/src/Persistent_cohomology/example/rips_persistence.cpp b/src/Persistent_cohomology/example/rips_persistence.cpp index c6378de7..d504798b 100644 --- a/src/Persistent_cohomology/example/rips_persistence.cpp +++ b/src/Persistent_cohomology/example/rips_persistence.cpp @@ -60,7 +60,7 @@ int main(int argc, char * argv[]) { program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); Points_off_reader off_reader(off_file_points); - Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Euclidean_distance()); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); // Construct the Rips complex in a Simplex Tree Simplex_tree simplex_tree; diff --git a/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp index b159c62e..554eeba6 100644 --- a/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp +++ b/src/Persistent_cohomology/example/rips_persistence_step_by_step.cpp @@ -82,7 +82,7 @@ int main(int argc, char * argv[]) { // Compute the proximity graph of the points Graph_t prox_graph = compute_proximity_graph(off_reader.get_point_cloud(), threshold - , Euclidean_distance()); + , Gudhi::Euclidean_distance()); // Construct the Rips complex in a Simplex Tree Simplex_tree st; diff --git a/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp index 63da9847..9618f278 100644 --- a/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp +++ b/src/Persistent_cohomology/example/rips_persistence_via_boundary_matrix.cpp @@ -70,7 +70,7 @@ int main(int argc, char * argv[]) { program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); Points_off_reader off_reader(off_file_points); - Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Euclidean_distance()); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); // Construct the Rips complex in a Simplex Tree Simplex_tree& st = *new Simplex_tree; diff --git a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp index 0ed3fddf..da418034 100644 --- a/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp +++ b/src/Persistent_cohomology/test/betti_numbers_unit_test.cpp @@ -9,8 +9,8 @@ #define BOOST_TEST_MODULE "betti_numbers" #include <boost/test/unit_test.hpp> -#include "gudhi/Simplex_tree.h" -#include "gudhi/Persistent_cohomology.h" +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> struct MiniSTOptions : Gudhi::Simplex_tree_options_full_featured { // Implicitly use 0 as filtration value for all simplices diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp index 6efd749e..e2e0bc71 100644 --- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp +++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test.cpp @@ -10,10 +10,10 @@ #define BOOST_TEST_MODULE "persistent_cohomology" #include <boost/test/unit_test.hpp> -#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/reader_utils.h" -#include "gudhi/Simplex_tree.h" -#include "gudhi/Persistent_cohomology.h" +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/reader_utils.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> using namespace Gudhi; using namespace Gudhi::persistent_cohomology; diff --git a/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp index 1a6e3296..6a7ecbec 100644 --- a/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp +++ b/src/Persistent_cohomology/test/persistent_cohomology_unit_test_multi_field.cpp @@ -9,11 +9,11 @@ #define BOOST_TEST_MODULE "persistent_cohomology_multi_field" #include <boost/test/unit_test.hpp> -#include "gudhi/graph_simplicial_complex.h" -#include "gudhi/reader_utils.h" -#include "gudhi/Simplex_tree.h" -#include "gudhi/Persistent_cohomology.h" -#include "gudhi/Persistent_cohomology/Multi_field.h" +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/reader_utils.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Persistent_cohomology/Multi_field.h> using namespace Gudhi; using namespace Gudhi::persistent_cohomology; diff --git a/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp b/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp index 3fd69ebc..a1db8910 100644 --- a/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp +++ b/src/Rips_complex/example/example_one_skeleton_rips_from_points.cpp @@ -27,7 +27,7 @@ int main() { // Init of a Rips complex from points // ---------------------------------------------------------------------------- double threshold = 12.0; - Rips_complex rips_complex_from_points(points, threshold, Euclidean_distance()); + Rips_complex rips_complex_from_points(points, threshold, Gudhi::Euclidean_distance()); Simplex_tree stree; rips_complex_from_points.create_complex(stree, 1); diff --git a/src/Rips_complex/example/example_rips_complex_from_off_file.cpp b/src/Rips_complex/example/example_rips_complex_from_off_file.cpp index a1e4e255..de2e4ea4 100644 --- a/src/Rips_complex/example/example_rips_complex_from_off_file.cpp +++ b/src/Rips_complex/example/example_rips_complex_from_off_file.cpp @@ -32,7 +32,7 @@ int main(int argc, char **argv) { // Init of a Rips complex from an OFF file // ---------------------------------------------------------------------------- Gudhi::Points_off_reader<Point> off_reader(off_file_name); - Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Euclidean_distance()); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); std::streambuf* streambufffer; std::ofstream ouput_file_stream; diff --git a/src/Rips_complex/test/test_rips_complex.cpp b/src/Rips_complex/test/test_rips_complex.cpp index ae68ba0d..fc2179f2 100644 --- a/src/Rips_complex/test/test_rips_complex.cpp +++ b/src/Rips_complex/test/test_rips_complex.cpp @@ -60,7 +60,7 @@ BOOST_AUTO_TEST_CASE(RIPS_DOC_OFF_file) { rips_threshold << "==========" << std::endl; Gudhi::Points_off_reader<Point> off_reader(off_file_name); - Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), rips_threshold, Euclidean_distance()); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), rips_threshold, Gudhi::Euclidean_distance()); const int DIMENSION_1 = 1; Simplex_tree st; @@ -89,10 +89,10 @@ BOOST_AUTO_TEST_CASE(RIPS_DOC_OFF_file) { std::cout << vertex << ","; vp.push_back(off_reader.get_point_cloud().at(vertex)); } - std::cout << ") - distance =" << Euclidean_distance()(vp.at(0), vp.at(1)) << + std::cout << ") - distance =" << Gudhi::Euclidean_distance()(vp.at(0), vp.at(1)) << " - filtration =" << st.filtration(f_simplex) << std::endl; BOOST_CHECK(vp.size() == 2); - BOOST_CHECK(are_almost_the_same(st.filtration(f_simplex), Euclidean_distance()(vp.at(0), vp.at(1)))); + BOOST_CHECK(are_almost_the_same(st.filtration(f_simplex), Gudhi::Euclidean_distance()(vp.at(0), vp.at(1)))); } } @@ -341,7 +341,7 @@ BOOST_AUTO_TEST_CASE(Rips_create_complex_throw) { rips_threshold << "==========" << std::endl; Gudhi::Points_off_reader<Point> off_reader(off_file_name); - Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), rips_threshold, Euclidean_distance()); + Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), rips_threshold, Gudhi::Euclidean_distance()); Simplex_tree stree; std::vector<int> simplex = {0, 1, 2}; diff --git a/src/Subsampling/example/example_choose_n_farthest_points.cpp b/src/Subsampling/example/example_choose_n_farthest_points.cpp index 533aba74..ebf631fc 100644 --- a/src/Subsampling/example/example_choose_n_farthest_points.cpp +++ b/src/Subsampling/example/example_choose_n_farthest_points.cpp @@ -19,7 +19,9 @@ int main(void) { K k; std::vector<Point_d> results; - Gudhi::subsampling::choose_n_farthest_points(k, points, 100, std::back_inserter(results)); + Gudhi::subsampling::choose_n_farthest_points(k, points, 100, + Gudhi::subsampling::random_starting_point, + std::back_inserter(results)); std::cout << "Before sparsification: " << points.size() << " points.\n"; std::cout << "After sparsification: " << results.size() << " points.\n"; diff --git a/src/Subsampling/example/example_custom_kernel.cpp b/src/Subsampling/example/example_custom_kernel.cpp index 25b5bf6c..2d42bdde 100644 --- a/src/Subsampling/example/example_custom_kernel.cpp +++ b/src/Subsampling/example/example_custom_kernel.cpp @@ -54,7 +54,9 @@ int main(void) { std::vector<Point_d> points = {0, 1, 2, 3}; std::vector<Point_d> results; - Gudhi::subsampling::choose_n_farthest_points(k, points, 2, std::back_inserter(results)); + Gudhi::subsampling::choose_n_farthest_points(k, points, 2, + Gudhi::subsampling::random_starting_point, + std::back_inserter(results)); std::cout << "Before sparsification: " << points.size() << " points.\n"; std::cout << "After sparsification: " << results.size() << " points.\n"; std::cout << "Result table: {" << results[0] << "," << results[1] << "}\n"; diff --git a/src/Subsampling/include/gudhi/choose_n_farthest_points.h b/src/Subsampling/include/gudhi/choose_n_farthest_points.h index 5e908090..86500b28 100644 --- a/src/Subsampling/include/gudhi/choose_n_farthest_points.h +++ b/src/Subsampling/include/gudhi/choose_n_farthest_points.h @@ -25,16 +25,9 @@ #include <boost/range.hpp> -#include <gudhi/Kd_tree_search.h> - -#include <gudhi/Clock.h> - -#include <CGAL/Search_traits.h> -#include <CGAL/Search_traits_adapter.h> -#include <CGAL/Fuzzy_sphere.h> +#include <gudhi/Null_output_iterator.h> #include <iterator> -#include <algorithm> // for sort #include <vector> #include <random> #include <limits> // for numeric_limits<> @@ -43,36 +36,51 @@ namespace Gudhi { namespace subsampling { +/** + * \ingroup subsampling + */ +enum : std::size_t { +/** + * Argument for `choose_n_farthest_points` to indicate that the starting point should be picked randomly. + */ + random_starting_point = std::size_t(-1) +}; + /** * \ingroup subsampling * \brief Subsample by a greedy strategy of iteratively adding the farthest point from the * current chosen point set to the subsampling. - * The iteration starts with the landmark `starting point`. + * The iteration starts with the landmark `starting point` or, if `starting point==random_starting_point`, with a random landmark. * \tparam Kernel must provide a type Kernel::Squared_distance_d which is a model of the * concept <a target="_blank" - * href="http://doc.cgal.org/latest/Kernel_d/classKernel__d_1_1Squared__distance__d.html">Kernel_d::Squared_distance_d</a> - * concept. - * It must also contain a public member 'squared_distance_d_object' of this type. + * href="http://doc.cgal.org/latest/Kernel_d/classKernel__d_1_1Squared__distance__d.html">Kernel_d::Squared_distance_d</a> (despite the name, taken from CGAL, this can be any kind of metric or proximity measure). + * It must also contain a public member `squared_distance_d_object()` that returns an object of this type. * \tparam Point_range Range whose value type is Kernel::Point_d. It must provide random-access * via `operator[]` and the points should be stored contiguously in memory. - * \tparam OutputIterator Output iterator whose value type is Kernel::Point_d. - * \details It chooses `final_size` points from a random access range `input_pts` and - * outputs it in the output iterator `output_it`. + * \tparam PointOutputIterator Output iterator whose value type is Kernel::Point_d. + * \tparam DistanceOutputIterator Output iterator for distances. + * \details It chooses `final_size` points from a random access range + * `input_pts` and outputs them in the output iterator `output_it`. It also + * outputs the distance from each of those points to the set of previous + * points in `dist_it`. * @param[in] k A kernel object. * @param[in] input_pts Const reference to the input points. * @param[in] final_size The size of the subsample to compute. * @param[in] starting_point The seed in the farthest point algorithm. - * @param[out] output_it The output iterator. + * @param[out] output_it The output iterator for points. + * @param[out] dist_it The optional output iterator for distances. * */ template < typename Kernel, typename Point_range, -typename OutputIterator> +typename PointOutputIterator, +typename DistanceOutputIterator = Null_output_iterator> void choose_n_farthest_points(Kernel const &k, Point_range const &input_pts, std::size_t final_size, std::size_t starting_point, - OutputIterator output_it) { + PointOutputIterator output_it, + DistanceOutputIterator dist_it = {}) { std::size_t nb_points = boost::size(input_pts); if (final_size > nb_points) final_size = nb_points; @@ -81,6 +89,14 @@ void choose_n_farthest_points(Kernel const &k, if (final_size < 1) return; + if (starting_point == random_starting_point) { + // Choose randomly the first landmark + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_int_distribution<std::size_t> dis(0, (input_pts.size() - 1)); + starting_point = dis(gen); + } + typename Kernel::Squared_distance_d sqdist = k.squared_distance_d_object(); std::size_t current_number_of_landmarks = 0; // counter for landmarks @@ -92,6 +108,7 @@ void choose_n_farthest_points(Kernel const &k, for (current_number_of_landmarks = 0; current_number_of_landmarks != final_size; current_number_of_landmarks++) { // curr_max_w at this point is the next landmark *output_it++ = input_pts[curr_max_w]; + *dist_it++ = dist_to_L[curr_max_w]; std::size_t i = 0; for (auto& p : input_pts) { double curr_dist = sqdist(p, *(std::begin(input_pts) + curr_max_w)); @@ -109,47 +126,6 @@ void choose_n_farthest_points(Kernel const &k, } } -/** - * \ingroup subsampling - * \brief Subsample by a greedy strategy of iteratively adding the farthest point from the - * current chosen point set to the subsampling. - * The iteration starts with a random landmark. - * \tparam Kernel must provide a type Kernel::Squared_distance_d which is a model of the - * concept <a target="_blank" - * href="http://doc.cgal.org/latest/Kernel_d/classKernel__d_1_1Squared__distance__d.html">Kernel_d::Squared_distance_d</a> - * concept. - * It must also contain a public member 'squared_distance_d_object' of this type. - * \tparam Point_range Range whose value type is Kernel::Point_d. It must provide random-access - * via `operator[]` and the points should be stored contiguously in memory. - * \tparam OutputIterator Output iterator whose value type is Kernel::Point_d. - * \details It chooses `final_size` points from a random access range `input_pts` and - * outputs it in the output iterator `output_it`. - * @param[in] k A kernel object. - * @param[in] input_pts Const reference to the input points. - * @param[in] final_size The size of the subsample to compute. - * @param[out] output_it The output iterator. - * - */ -template < typename Kernel, -typename Point_range, -typename OutputIterator> -void choose_n_farthest_points(Kernel const& k, - Point_range const &input_pts, - unsigned final_size, - OutputIterator output_it) { - // Tests to the limit - if ((final_size < 1) || (input_pts.size() == 0)) - return; - - // Choose randomly the first landmark - std::random_device rd; - std::mt19937 gen(rd()); - std::uniform_int_distribution<> dis(0, (input_pts.size() - 1)); - std::size_t starting_point = dis(gen); - - choose_n_farthest_points(k, input_pts, final_size, starting_point, output_it); -} - } // namespace subsampling } // namespace Gudhi diff --git a/src/Subsampling/test/test_choose_n_farthest_points.cpp b/src/Subsampling/test/test_choose_n_farthest_points.cpp index 0bc0dff4..ee9d4c77 100644 --- a/src/Subsampling/test/test_choose_n_farthest_points.cpp +++ b/src/Subsampling/test/test_choose_n_farthest_points.cpp @@ -25,7 +25,7 @@ // #endif #define BOOST_TEST_DYN_LINK -#define BOOST_TEST_MODULE "witness_complex_points" +#define BOOST_TEST_MODULE Subsampling - test choose_n_farthest_points #include <boost/test/unit_test.hpp> #include <boost/mpl/list.hpp> @@ -56,7 +56,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_choose_farthest_point, Kernel, list_of_tested landmarks.clear(); Kernel k; - Gudhi::subsampling::choose_n_farthest_points(k, points, 100, std::back_inserter(landmarks)); + Gudhi::subsampling::choose_n_farthest_points(k, points, 100, Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); BOOST_CHECK(landmarks.size() == 100); for (auto landmark : landmarks) @@ -70,34 +70,45 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_choose_farthest_point_limits, Kernel, list_of typedef typename Kernel::FT FT; typedef typename Kernel::Point_d Point_d; std::vector< Point_d > points, landmarks; + std::vector< FT > distances; landmarks.clear(); Kernel k; // Choose -1 farthest points in an empty point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, -1, std::back_inserter(landmarks)); + Gudhi::subsampling::choose_n_farthest_points(k, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); BOOST_CHECK(landmarks.size() == 0); - landmarks.clear(); + landmarks.clear(); distances.clear(); // Choose 0 farthest points in an empty point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, 0, std::back_inserter(landmarks)); + Gudhi::subsampling::choose_n_farthest_points(k, points, 0, -1, std::back_inserter(landmarks), std::back_inserter(distances)); BOOST_CHECK(landmarks.size() == 0); - landmarks.clear(); + landmarks.clear(); distances.clear(); // Choose 1 farthest points in an empty point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, 1, std::back_inserter(landmarks)); + Gudhi::subsampling::choose_n_farthest_points(k, points, 1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); BOOST_CHECK(landmarks.size() == 0); - landmarks.clear(); + landmarks.clear(); distances.clear(); std::vector<FT> point({0.0, 0.0, 0.0, 0.0}); points.push_back(Point_d(point.begin(), point.end())); - // Choose -1 farthest points in an empty point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, -1, std::back_inserter(landmarks)); - BOOST_CHECK(landmarks.size() == 1); - landmarks.clear(); + // Choose -1 farthest points in a one point cloud + Gudhi::subsampling::choose_n_farthest_points(k, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + BOOST_CHECK(landmarks.size() == 1 && distances.size() == 1); + BOOST_CHECK(distances[0] == std::numeric_limits<FT>::infinity()); + landmarks.clear(); distances.clear(); // Choose 0 farthest points in a one point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, 0, std::back_inserter(landmarks)); - BOOST_CHECK(landmarks.size() == 0); - landmarks.clear(); + Gudhi::subsampling::choose_n_farthest_points(k, points, 0, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + BOOST_CHECK(landmarks.size() == 0 && distances.size() == 0); + landmarks.clear(); distances.clear(); // Choose 1 farthest points in a one point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, 1, std::back_inserter(landmarks)); - BOOST_CHECK(landmarks.size() == 1); - landmarks.clear(); + Gudhi::subsampling::choose_n_farthest_points(k, points, 1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + BOOST_CHECK(landmarks.size() == 1 && distances.size() == 1); + BOOST_CHECK(distances[0] == std::numeric_limits<FT>::infinity()); + landmarks.clear(); distances.clear(); + std::vector<FT> point2({1.0, 0.0, 0.0, 0.0}); + points.push_back(Point_d(point2.begin(), point2.end())); + // Choose all farthest points in a one point cloud + Gudhi::subsampling::choose_n_farthest_points(k, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + BOOST_CHECK(landmarks.size() == 2 && distances.size() == 2); + BOOST_CHECK(distances[0] == std::numeric_limits<FT>::infinity()); + BOOST_CHECK(distances[1] == 1); + landmarks.clear(); distances.clear(); } diff --git a/src/common/include/gudhi/Null_output_iterator.h b/src/common/include/gudhi/Null_output_iterator.h new file mode 100644 index 00000000..42e6e449 --- /dev/null +++ b/src/common/include/gudhi/Null_output_iterator.h @@ -0,0 +1,48 @@ +/* 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): Marc Glisse + * + * Copyright (C) 2017 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef NULL_OUTPUT_ITERATOR_H_ +#define NULL_OUTPUT_ITERATOR_H_ + +#include <iterator> + +namespace Gudhi { + +/** An output iterator that ignores whatever it is given. */ +struct Null_output_iterator { + typedef std::output_iterator_tag iterator_category; + typedef void value_type; + typedef void difference_type; + typedef void pointer; + typedef void reference; + + Null_output_iterator& operator++() {return *this;} + Null_output_iterator operator++(int) {return *this;} + struct proxy { + template<class T> + proxy& operator=(T&&){return *this;} + }; + proxy operator*()const{return {};} +}; +} // namespace Gudhi + +#endif // NULL_OUTPUT_ITERATOR_H_ diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index 22747637..f6e2ab5a 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -27,6 +27,8 @@ #include <type_traits> // for std::decay #include <iterator> // for std::begin, std::end +namespace Gudhi { + /** @file * @brief Global distance functions */ @@ -48,4 +50,6 @@ class Euclidean_distance { } }; +} // namespace Gudhi + #endif // DISTANCE_FUNCTIONS_H_ diff --git a/src/cython/include/Rips_complex_interface.h b/src/cython/include/Rips_complex_interface.h index 1879bd74..6d813f4a 100644 --- a/src/cython/include/Rips_complex_interface.h +++ b/src/cython/include/Rips_complex_interface.h @@ -49,7 +49,7 @@ class Rips_complex_interface { if (euclidean) { // Rips construction where values is a vector of points rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(values, threshold, - Euclidean_distance()); + Gudhi::Euclidean_distance()); } else { // Rips construction where values is a distance matrix rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(values, threshold); @@ -61,7 +61,8 @@ class Rips_complex_interface { // Rips construction where file_name is an OFF file Gudhi::Points_off_reader<Point_d> off_reader(file_name); rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(off_reader.get_point_cloud(), - threshold, Euclidean_distance()); + threshold, + Gudhi::Euclidean_distance()); } else { // Rips construction where values is a distance matrix Distance_matrix distances = diff --git a/src/cython/include/Subsampling_interface.h b/src/cython/include/Subsampling_interface.h index 1c6032c0..b0f4a50a 100644 --- a/src/cython/include/Subsampling_interface.h +++ b/src/cython/include/Subsampling_interface.h @@ -46,7 +46,7 @@ std::vector<std::vector<double>> subsampling_n_farthest_points(const std::vector unsigned nb_points) { std::vector<std::vector<double>> landmarks; Subsampling_dynamic_kernel k; - choose_n_farthest_points(k, points, nb_points, std::back_inserter(landmarks)); + choose_n_farthest_points(k, points, nb_points, random_starting_point, std::back_inserter(landmarks)); return landmarks; } |