diff options
Diffstat (limited to 'src/Rips_complex')
-rw-r--r-- | src/Rips_complex/doc/Intro_rips_complex.h | 24 | ||||
-rw-r--r-- | src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 115 |
2 files changed, 99 insertions, 40 deletions
diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index b2840686..3888ec8f 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -64,7 +64,7 @@ namespace rips_complex { * And so on for simplex (0,1,2,3). * * If the Rips_complex interfaces are not detailed enough for your need, please refer to - * <a href="_persistent_cohomology_2rips_persistence_step_by_step_8cpp-example.html"> + * <a href="rips_persistence_step_by_step_8cpp-example.html"> * rips_persistence_step_by_step.cpp</a> example, where the constructions of the graph and * the Simplex_tree are more detailed. * @@ -111,7 +111,7 @@ namespace rips_complex { * * Then, it is asked to display information about the simplicial complex. * - * \include Rips_complex/example_one_skeleton_rips_from_points.cpp + * \include example_one_skeleton_rips_from_points.cpp * * When launching (Rips maximal distance between 2 points is 12.0, is expanded * until dimension 1 - one skeleton graph in other words): @@ -121,7 +121,7 @@ namespace rips_complex { * * the program output is: * - * \include Rips_complex/one_skeleton_rips_for_doc.txt + * \include one_skeleton_rips_for_doc.txt * * \subsection ripsoffexample Example from OFF file * @@ -132,7 +132,7 @@ namespace rips_complex { * * Then, it is asked to display information about the Rips complex. * - * \include Rips_complex/example_rips_complex_from_off_file.cpp + * \include example_rips_complex_from_off_file.cpp * * When launching: * @@ -141,7 +141,7 @@ namespace rips_complex { * * the program output is: * - * \include Rips_complex/full_skeleton_rips_for_doc.txt + * \include full_skeleton_rips_for_doc.txt * * * \subsection sparseripspointscloudexample Example of a sparse Rips from a point cloud @@ -149,7 +149,7 @@ namespace rips_complex { * This example builds the full sparse Rips of a set of 2D Euclidean points, then prints some minimal * information about the complex. * - * \include Rips_complex/example_sparse_rips.cpp + * \include example_sparse_rips.cpp * * When launching: * @@ -172,7 +172,7 @@ namespace rips_complex { * * Then, it is asked to display information about the simplicial complex. * - * \include Rips_complex/example_one_skeleton_rips_from_distance_matrix.cpp + * \include example_one_skeleton_rips_from_distance_matrix.cpp * * When launching (Rips maximal distance between 2 points is 1.0, is expanded until dimension 1 - one skeleton graph * with other words): @@ -182,7 +182,7 @@ namespace rips_complex { * * the program output is: * - * \include Rips_complex/one_skeleton_rips_for_doc.txt + * \include one_skeleton_rips_for_doc.txt * * \subsection ripscsvdistanceexample Example from a distance matrix read in a csv file * @@ -192,7 +192,7 @@ namespace rips_complex { * * Then, it is asked to display information about the Rips complex. * - * \include Rips_complex/example_rips_complex_from_csv_distance_matrix_file.cpp + * \include example_rips_complex_from_csv_distance_matrix_file.cpp * * When launching: * @@ -201,7 +201,7 @@ namespace rips_complex { * * the program output is: * - * \include Rips_complex/full_skeleton_rips_for_doc.txt + * \include full_skeleton_rips_for_doc.txt * * * \section ripscorrelationematrix Correlation matrix @@ -213,7 +213,7 @@ namespace rips_complex { * * Then, it is asked to display information about the simplicial complex. * - * \include Rips_complex/example_one_skeleton_rips_from_correlation_matrix.cpp + * \include example_one_skeleton_rips_from_correlation_matrix.cpp * * When launching: * @@ -222,7 +222,7 @@ namespace rips_complex { * * the program output is: * - * \include Rips_complex/one_skeleton_rips_from_correlation_matrix_for_doc.txt + * \include one_skeleton_rips_from_correlation_matrix_for_doc.txt * * All the other constructions discussed for Rips complex for distance matrix can be also performed for Rips complexes * construction from correlation matrices. diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index a5501004..7ae7b317 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -15,12 +15,71 @@ #include <gudhi/graph_simplicial_complex.h> #include <gudhi/choose_n_farthest_points.h> -#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/graph_traits.hpp> #include <boost/range/metafunctions.hpp> +#include <boost/iterator/counting_iterator.hpp> #include <vector> namespace Gudhi { +namespace rips_complex { +// A custom graph class, because boost::adjacency_list does not conveniently allow to choose vertex descriptors +template <class Vertex_handle, class Filtration_value> +struct Graph { + typedef std::vector<Vertex_handle> VList; + typedef std::vector<std::tuple<Vertex_handle, Vertex_handle, Filtration_value>> EList; + typedef typename VList::const_iterator vertex_iterator; + typedef boost::counting_iterator<std::size_t> edge_iterator; + VList vlist; + EList elist; +}; +template <class Vertex_handle, class Filtration_value> +void add_vertex(Vertex_handle v, Graph<Vertex_handle, Filtration_value>&g) { g.vlist.push_back(v); } +template <class Vertex_handle, class Filtration_value> +void add_edge(Vertex_handle u, Vertex_handle v, Filtration_value f, Graph<Vertex_handle, Filtration_value>&g) { g.elist.emplace_back(u, v, f); } +template <class Vertex_handle, class Filtration_value> +std::size_t num_vertices(Graph<Vertex_handle, Filtration_value> const&g) { return g.vlist.size(); } +template <class Vertex_handle, class Filtration_value> +std::size_t num_edges(Graph<Vertex_handle, Filtration_value> const&g) { return g.elist.size(); } +template <class Vertex_handle, class Filtration_value, class Iter = typename Graph<Vertex_handle, Filtration_value>::vertex_iterator> +std::pair<Iter, Iter> +vertices(Graph<Vertex_handle, Filtration_value> const&g) { + return { g.vlist.begin(), g.vlist.end() }; +} +template <class Vertex_handle, class Filtration_value> +std::pair<boost::counting_iterator<std::size_t>, boost::counting_iterator<std::size_t>> +edges(Graph<Vertex_handle, Filtration_value> const&g) { + typedef boost::counting_iterator<std::size_t> I; + return { I(0), I(g.elist.size()) }; +} +template <class Vertex_handle, class Filtration_value> +Vertex_handle source(std::size_t e, Graph<Vertex_handle, Filtration_value> const&g) { return std::get<0>(g.elist[e]); } +template <class Vertex_handle, class Filtration_value> +Vertex_handle target(std::size_t e, Graph<Vertex_handle, Filtration_value> const&g) { return std::get<1>(g.elist[e]); } +template <class Vertex_handle, class Filtration_value> +Filtration_value get(vertex_filtration_t, Graph<Vertex_handle, Filtration_value> const&, Vertex_handle) { return 0; } +template <class Vertex_handle, class Filtration_value> +Filtration_value get(edge_filtration_t, Graph<Vertex_handle, Filtration_value> const&g, std::size_t e) { return std::get<2>(g.elist[e]); } +} // namespace rips_complex +} // namespace Gudhi +namespace boost { +template <class Vertex_handle, class Filtration_value> +struct graph_traits<Gudhi::rips_complex::Graph<Vertex_handle, Filtration_value>> { + typedef Gudhi::rips_complex::Graph<Vertex_handle, Filtration_value> G; + struct traversal_category : vertex_list_graph_tag, edge_list_graph_tag {}; + typedef Vertex_handle vertex_descriptor; + typedef typename G::vertex_iterator vertex_iterator; + typedef std::size_t vertices_size_type; + typedef std::size_t edge_descriptor; + typedef typename G::edge_iterator edge_iterator; + typedef std::size_t edges_size_type; + 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 for instance. +} + +namespace Gudhi { namespace rips_complex { @@ -45,12 +104,8 @@ template <typename Filtration_value> class Sparse_rips_complex { private: // TODO(MG): use a different graph where we know we can safely insert in parallel. - typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS, - boost::property<vertex_filtration_t, Filtration_value>, - boost::property<edge_filtration_t, Filtration_value>> - Graph; - typedef int Vertex_handle; + typedef rips_complex::Graph<Vertex_handle, Filtration_value> Graph; public: /** \brief Sparse_rips_complex constructor from a list of points. @@ -63,10 +118,11 @@ class Sparse_rips_complex { * */ template <typename RandomAccessPointRange, typename Distance> - Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon, Filtration_value mini=-std::numeric_limits<Filtration_value>::infinity(), Filtration_value maxi=std::numeric_limits<Filtration_value>::infinity()) + Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double const epsilon, Filtration_value const mini=-std::numeric_limits<Filtration_value>::infinity(), Filtration_value const maxi=std::numeric_limits<Filtration_value>::infinity()) : 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 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); @@ -83,7 +139,7 @@ class Sparse_rips_complex { * @param[in] maxi Maximal filtration value. Ignore anything above this scale. */ template <typename DistanceMatrix> - Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon, Filtration_value mini=-std::numeric_limits<Filtration_value>::infinity(), Filtration_value maxi=std::numeric_limits<Filtration_value>::infinity()) + Sparse_rips_complex(const DistanceMatrix& distance_matrix, double const epsilon, Filtration_value const mini=-std::numeric_limits<Filtration_value>::infinity(), Filtration_value const maxi=std::numeric_limits<Filtration_value>::infinity()) : Sparse_rips_complex(boost::irange<Vertex_handle>(0, boost::size(distance_matrix)), [&](Vertex_handle i, Vertex_handle j) { return (i==j) ? 0 : (i<j) ? distance_matrix[j][i] : distance_matrix[i][j]; }, epsilon, mini, maxi) {} @@ -99,7 +155,7 @@ class Sparse_rips_complex { * */ template <typename SimplicialComplexForRips> - void create_complex(SimplicialComplexForRips& complex, int dim_max) { + void create_complex(SimplicialComplexForRips& complex, int const dim_max) { GUDHI_CHECK(complex.num_vertices() == 0, std::invalid_argument("Sparse_rips_complex::create_complex - simplicial complex is not empty")); @@ -108,17 +164,17 @@ class Sparse_rips_complex { complex.expansion(dim_max); return; } - const int n = boost::size(params); - std::vector<Filtration_value> lambda(n); + const std::size_t 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(std::size_t 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){ auto filt = complex.filtration(sh); - auto mini = filt * cst; + auto min_f = filt * cst; for(auto v : complex.simplex_vertex_range(sh)){ - if(lambda[v] < mini) + if(lambda[v] < min_f) return true; // v died before this simplex could be born } return false; @@ -129,32 +185,34 @@ class Sparse_rips_complex { private: // PointRange must be random access. template <typename Distance> - void compute_sparse_graph(Distance& dist, double epsilon, Filtration_value mini, Filtration_value maxi) { + void compute_sparse_graph(Distance& dist, double const epsilon, Filtration_value const mini, Filtration_value const maxi) { const auto& points = sorted_points; // convenience alias - const int n = boost::size(points); + std::size_t n = boost::size(points); double cst = epsilon * (1 - epsilon) / 2; - graph_.~Graph(); - new (&graph_) Graph(n); - // for(auto v : vertices(g)) // doesn't work :-( - typename boost::graph_traits<Graph>::vertex_iterator v_i, v_e; - for (std::tie(v_i, v_e) = vertices(graph_); v_i != v_e; ++v_i) { - auto v = *v_i; - // This whole loop might not be necessary, leave it until someone investigates if it is safe to remove. - put(vertex_filtration_t(), graph_, v, 0); + max_v = -1; // Useful for the size of the map lambda. + for (std::size_t 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 (std::size_t i = 0; i < n; ++i) { auto&& pi = points[i]; auto li = params[i]; - if (li < mini) break; - 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 (std::size_t j = i + 1; j < n; ++j) { auto&& pj = points[j]; auto d = dist(pi, pj); auto lj = params[j]; - if (lj < mini) break; GUDHI_CHECK(lj <= li, "Bad furthest point sorting"); Filtration_value alpha; @@ -178,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; |