diff options
Diffstat (limited to 'src/Rips_complex/include/gudhi/Sparse_rips_complex.h')
| -rw-r--r-- | src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 236 |
1 files changed, 119 insertions, 117 deletions
diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 503a783a..1a9d6ebb 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -32,138 +32,140 @@ #include <vector> - namespace Gudhi { namespace rips_complex { -// The whole interface is copied on Rips_complex. A redesign should be discussed with all complex creation classes in mind. +// The whole interface is copied on Rips_complex. A redesign should be discussed with all complex creation classes in +// mind. /** * \class Sparse_rips_complex * \brief Sparse Rips complex data structure. - * + * * \ingroup rips_complex - * + * * \details - * This class is used to construct a sparse \f$(1+O(\epsilon))\f$-approximation of `Rips_complex`, i.e. a filtered simplicial complex that is multiplicatively \f$(1+O(\epsilon))\f$-interleaved with the Rips filtration. - * + * This class is used to construct a sparse \f$(1+O(\epsilon))\f$-approximation of `Rips_complex`, i.e. a filtered + * simplicial complex that is multiplicatively \f$(1+O(\epsilon))\f$-interleaved with the Rips filtration. + * * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex. */ -template<typename Filtration_value> +template <typename Filtration_value> class Sparse_rips_complex { - private: - // TODO: use a different graph where we know we can safely insert in parallel. - typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS - , boost::property<vertex_filtration_t, Filtration_value> - , boost::property<edge_filtration_t, Filtration_value>> Graph; - - typedef int Vertex_handle; - - public: - /** \brief Sparse_rips_complex constructor from a list of points. - * - * @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. - * - */ - template<typename RandomAccessPointRange, typename Distance > - Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) { - GUDHI_CHECK(epsilon > 0, "epsilon must be positive"); - std::vector<Vertex_handle> sorted_points; - std::vector<Filtration_value> params; - auto dist_fun = [&](Vertex_handle i, Vertex_handle j){return distance(points[i], points[j]);}; - Ker<decltype(dist_fun)> kernel(dist_fun); - subsampling::choose_n_farthest_points(kernel, boost::irange<Vertex_handle>(0, boost::size(points)), -1, -1, std::back_inserter(sorted_points), std::back_inserter(params)); - compute_sparse_graph(sorted_points, params, dist_fun, epsilon); - } - - /** \brief Sparse_rips_complex constructor from a distance matrix. - * - * @param[in] distance_matrix Range of range of distances. - * `distance_matrix[i][j]` returns the distance between points \f$i\f$ and - * \f$j\f$ as long as \f$ 0 \leqslant i < j \leqslant - * distance\_matrix.size().\f$ - * @param[in] epsilon Approximation parameter. epsilon must be positive. - */ - template<typename DistanceMatrix> - Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon) - : Sparse_rips_complex( - boost::irange<Vertex_handle>(0, boost::size(distance_matrix)), - [&](Vertex_handle i, Vertex_handle j){return distance_matrix[j][i];}, - epsilon) {} - - /** \brief Fills the simplicial complex with the sparse Rips graph and - * expands it with all the cliques, stopping at a given maximal dimension. - * - * \tparam SimplicialComplexForRips must meet `SimplicialComplexForRips` concept. - * - * @param[in] complex the complex to fill - * @param[in] dim_max maximal dimension of the simplicial complex. - * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. - * - */ - template <typename SimplicialComplexForRips> - void create_complex(SimplicialComplexForRips& complex, int dim_max) { - GUDHI_CHECK(complex.num_vertices() == 0, - std::invalid_argument("Sparse_rips_complex::create_complex - simplicial complex is not empty")); - - complex.insert_graph(graph_); - complex.expansion(dim_max); - } - - private: - // choose_n_farthest_points wants the distance function in this form... - template <class Distance> - struct Ker { - typedef std::size_t Point_d; // index into point range - Ker(Distance& d) : dist (d) {} - // Despite the name, this is not squared... - typedef Distance Squared_distance_d; - Squared_distance_d& squared_distance_d_object() const { return dist; } - Distance& dist; - }; - - // PointRange must be random access. - template<typename PointRange, typename ParamRange, typename Distance> - void compute_sparse_graph(const PointRange& points, const ParamRange& params, Distance& dist, double epsilon) { - const int n = boost::size(points); - 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); - } - - // TODO: - // - make it parallel - // - only test near-enough neighbors - for(int i = 0; i < n; ++i) - for(int j = i + 1; j < n; ++j){ - auto&& pi = points[i]; - auto&& pj = points[j]; - auto d = dist(pi, pj); - auto li = params[i]; - auto lj = params[j]; - GUDHI_CHECK(lj <= li, "Bad furthest point sorting"); - Filtration_value alpha; - - // The paper has d/2 and d-lj/e to match the Cech, but we use doubles to match the Rips - if(d * epsilon <= 2 * lj) - alpha = d; - else if(d * epsilon <= li + lj && (epsilon >= 1 || d * epsilon <= lj * (1 + 1 / (1 - epsilon)))) - alpha = (d - lj / epsilon) * 2; - else continue; - - add_edge(pi, pj, alpha, graph_); - } + private: + // TODO: use a different graph where we know we can safely insert in parallel. + typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS, + boost::property<vertex_filtration_t, Filtration_value>, + boost::property<edge_filtration_t, Filtration_value>> + Graph; + + typedef int Vertex_handle; + + public: + /** \brief Sparse_rips_complex constructor from a list of points. + * + * @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. + * + */ + template <typename RandomAccessPointRange, typename Distance> + Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) { + GUDHI_CHECK(epsilon > 0, "epsilon must be positive"); + std::vector<Vertex_handle> sorted_points; + std::vector<Filtration_value> params; + auto dist_fun = [&](Vertex_handle i, Vertex_handle j) { return distance(points[i], points[j]); }; + Ker<decltype(dist_fun)> kernel(dist_fun); + subsampling::choose_n_farthest_points(kernel, boost::irange<Vertex_handle>(0, boost::size(points)), -1, -1, + std::back_inserter(sorted_points), std::back_inserter(params)); + compute_sparse_graph(sorted_points, params, dist_fun, epsilon); + } + + /** \brief Sparse_rips_complex constructor from a distance matrix. + * + * @param[in] distance_matrix Range of range of distances. + * `distance_matrix[i][j]` returns the distance between points \f$i\f$ and + * \f$j\f$ as long as \f$ 0 \leqslant i < j \leqslant + * distance\_matrix.size().\f$ + * @param[in] epsilon Approximation parameter. epsilon must be positive. + */ + template <typename DistanceMatrix> + Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon) + : Sparse_rips_complex(boost::irange<Vertex_handle>(0, boost::size(distance_matrix)), + [&](Vertex_handle i, Vertex_handle j) { return distance_matrix[j][i]; }, epsilon) {} + + /** \brief Fills the simplicial complex with the sparse Rips graph and + * expands it with all the cliques, stopping at a given maximal dimension. + * + * \tparam SimplicialComplexForRips must meet `SimplicialComplexForRips` concept. + * + * @param[in] complex the complex to fill + * @param[in] dim_max maximal dimension of the simplicial complex. + * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. + * + */ + template <typename SimplicialComplexForRips> + void create_complex(SimplicialComplexForRips& complex, int dim_max) { + GUDHI_CHECK(complex.num_vertices() == 0, + std::invalid_argument("Sparse_rips_complex::create_complex - simplicial complex is not empty")); + + complex.insert_graph(graph_); + complex.expansion(dim_max); + } + + private: + // choose_n_farthest_points wants the distance function in this form... + template <class Distance> + struct Ker { + typedef std::size_t Point_d; // index into point range + Ker(Distance& d) : dist(d) {} + // Despite the name, this is not squared... + typedef Distance Squared_distance_d; + Squared_distance_d& squared_distance_d_object() const { return dist; } + Distance& dist; + }; + + // PointRange must be random access. + template <typename PointRange, typename ParamRange, typename Distance> + void compute_sparse_graph(const PointRange& points, const ParamRange& params, Distance& dist, double epsilon) { + const int n = boost::size(points); + 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); + } + + // TODO: + // - make it parallel + // - only test near-enough neighbors + for (int i = 0; i < n; ++i) + for (int j = i + 1; j < n; ++j) { + auto&& pi = points[i]; + auto&& pj = points[j]; + auto d = dist(pi, pj); + auto li = params[i]; + auto lj = params[j]; + GUDHI_CHECK(lj <= li, "Bad furthest point sorting"); + Filtration_value alpha; + + // The paper has d/2 and d-lj/e to match the Cech, but we use doubles to match the Rips + if (d * epsilon <= 2 * lj) + alpha = d; + else if (d * epsilon <= li + lj && (epsilon >= 1 || d * epsilon <= lj * (1 + 1 / (1 - epsilon)))) + alpha = (d - lj / epsilon) * 2; + else + continue; + + add_edge(pi, pj, alpha, graph_); } + } - Graph graph_; + Graph graph_; }; } // namespace rips_complex |