diff options
author | Marc Glisse <marc.glisse@inria.fr> | 2021-04-17 17:31:23 +0200 |
---|---|---|
committer | Marc Glisse <marc.glisse@inria.fr> | 2021-04-17 18:52:37 +0200 |
commit | 71337179d95d1e330902b431907cb07698abcdc9 (patch) | |
tree | 542ef0cd93a40edf2c9414460855f4e0336e4a57 /src | |
parent | 20bee15a2e7dc68deb3141ebab7a30a3edcfb401 (diff) |
Safely drop some vertices in sparse Rips
Diffstat (limited to 'src')
-rw-r--r-- | src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 48 | ||||
-rwxr-xr-x | src/python/test/test_rips_complex.py | 2 |
2 files changed, 23 insertions, 27 deletions
diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 30afb1d0..28031e68 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -76,7 +76,7 @@ struct graph_traits<Gudhi::rips_complex::Graph<Vertex_handle, Filtration_value>> typedef directed_tag directed_category; typedef disallow_parallel_edge_tag edge_parallel_category; }; -// Etc, since we don't expose this graph to the world, we know we are not going to query property_traits. +// Etc, since we don't expose this graph to the world, we know we are not going to query property_traits for instance. } namespace Gudhi { @@ -113,8 +113,7 @@ class Sparse_rips_complex { * @param[in] points Range of points. * @param[in] distance Distance function that returns a `Filtration_value` from 2 given points. * @param[in] epsilon Approximation parameter. epsilon must be positive. - * @param[in] mini Not implemented yet, and broken in previous versions. Minimal filtration value. - * Ignore anything below this scale. This is a less efficient version of `Gudhi::subsampling::sparsify_point_set()`. + * @param[in] mini Minimal filtration value. Ignore anything below this scale. This is a less efficient version of `Gudhi::subsampling::sparsify_point_set()`. * @param[in] maxi Maximal filtration value. Ignore anything above this scale. * */ @@ -123,7 +122,7 @@ class Sparse_rips_complex { : epsilon_(epsilon) { GUDHI_CHECK(epsilon > 0, "epsilon must be positive"); auto dist_fun = [&](Vertex_handle i, Vertex_handle j) { return distance(points[i], points[j]); }; - // TODO: stop choose_n_farthest_points once it reaches mini? Then the graph vertices would not be [0, ..., n-1] which complicates things. + // TODO: stop choose_n_farthest_points once it reaches mini or 0? subsampling::choose_n_farthest_points(dist_fun, boost::irange<Vertex_handle>(0, boost::size(points)), -1, -1, std::back_inserter(sorted_points), std::back_inserter(params)); compute_sparse_graph(dist_fun, epsilon, mini, maxi); @@ -165,10 +164,10 @@ class Sparse_rips_complex { complex.expansion(dim_max); return; } - const int n = boost::size(params); - std::vector<Filtration_value> lambda(n); + const Vertex_handle n = num_vertices(graph_); + std::vector<Filtration_value> lambda(max_v + 1); // lambda[original_order]=params[sorted_order] - for(int i=0;i<n;++i) + for(Vertex_handle i=0;i<n;++i) lambda[sorted_points[i]] = params[i]; double cst = epsilon_ * (1 - epsilon_) / 2; auto block = [cst,&complex,&lambda](typename SimplicialComplexForRips::Simplex_handle sh){ @@ -188,36 +187,32 @@ class Sparse_rips_complex { template <typename Distance> void compute_sparse_graph(Distance& dist, double epsilon, Filtration_value mini, Filtration_value maxi) { const auto& points = sorted_points; // convenience alias - const int n = boost::size(points); + Vertex_handle n = boost::size(points); double cst = epsilon * (1 - epsilon) / 2; - graph_.~Graph(); - new (&graph_) Graph(); - for (int i = 0; i < n; ++i) { - add_vertex(i, graph_); + max_v = -1; // Useful for the size of the map lambda. + for (Vertex_handle i = 0; i < n; ++i) { + if ((params[i] < mini || params[i] <= 0) && i != 0) break; + // The parameter of the first point is not very meaningful, it is supposed to be infinite, + // but if the type does not support it... + // It would be better to do this reduction of the number of points earlier, around choose_n_farthest_points. + add_vertex(points[i], graph_); + max_v = std::max(max_v, points[i]); } + n = num_vertices(graph_); // TODO(MG): // - make it parallel // - only test near-enough neighbors - for (int i = 0; i < n; ++i) { + for (Vertex_handle i = 0; i < n; ++i) { auto&& pi = points[i]; auto li = params[i]; - // FIXME: see below about mini. It might be ok to uncomment just this one, but it requires a proof. - // if ((li < mini || li <= 0) && i != 0) break; - if (li <= 0 && i != 0) break; - // The parameter of the first point is not very meaningful, it is supposed to be infinite, - // but if the type does not support it... - // Points with multiplicity get connected to their first representative, no need to handle - // the redundant ones in the outer loop. - for (int j = i + 1; j < n; ++j) { + // If we inserted all the points, points with multiplicity would get connected to their first representative, + // no need to handle the redundant ones in the outer loop. + // if (li <= 0 && i != 0) break; + for (Vertex_handle j = i + 1; j < n; ++j) { auto&& pj = points[j]; auto d = dist(pi, pj); auto lj = params[j]; - // FIXME: It would make sense to ignore the points with low param completely, but the current graph type we are - // using implicitly inserts all the vertices 0 ... n-1, so this would create isolated vertices, which is bad. - // If we do end up ignoring those points, we should do it early, around choose_n_farthest_points. But be careful - // that the size of lambda should reflect the original number of points then. - // if (lj < mini) break; GUDHI_CHECK(lj <= li, "Bad furthest point sorting"); Filtration_value alpha; @@ -241,6 +236,7 @@ class Sparse_rips_complex { Graph graph_; double epsilon_; + Vertex_handle max_v; // Because of the arbitrary split between constructor and create_complex // sorted_points[sorted_order]=original_order std::vector<Vertex_handle> sorted_points; diff --git a/src/python/test/test_rips_complex.py b/src/python/test/test_rips_complex.py index cae21435..a2f43a1b 100755 --- a/src/python/test/test_rips_complex.py +++ b/src/python/test/test_rips_complex.py @@ -152,5 +152,5 @@ def test_sparse_with_multiplicity(): ] rips = RipsComplex(points=points, sparse=0.01) simplex_tree = rips.create_simplex_tree(max_dimension=2) - assert simplex_tree.num_simplices() == 25 + assert simplex_tree.num_simplices() == 7 diag = simplex_tree.persistence() |