From fe318b14148c20541d76649a1de76bdf6f7eaeab Mon Sep 17 00:00:00 2001 From: glisse Date: Tue, 16 Jan 2018 16:03:29 +0000 Subject: Initial commit of sparse Rips, based on Rips_complex. git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/sparserips-glisse@3137 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: f0550a0a24a25e8d9ab16535eb4364e8dc4228e6 --- .../include/gudhi/Sparse_rips_complex.h | 169 +++++++++++++++++++++ 1 file changed, 169 insertions(+) create mode 100644 src/Rips_complex/include/gudhi/Sparse_rips_complex.h (limited to 'src/Rips_complex/include') diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h new file mode 100644 index 00000000..c5378b6e --- /dev/null +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -0,0 +1,169 @@ +/* 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) 2018 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 . + */ + +#ifndef SPARSE_RIPS_COMPLEX_H_ +#define SPARSE_RIPS_COMPLEX_H_ + +#include +#include +#include + +#include +#include + +#include + + +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. + +/** + * \class Sparse_rips_complex + * \brief Sparse Rips complex data structure. + * + * \ingroup rips_complex + * + * \details + * This class is used to construct a sparse \f$(1+\epsilon)\f$-approximation of `Rips_complex`. + * + * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex. + */ +template +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::property> 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. + * + */ + template + Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) { + std::vector sorted_points; + std::vector params; + auto dist_fun = [&](Vertex_handle i, Vertex_handle j){return distance(points[i], points[j]);}; + Ker kernel(dist_fun); + subsampling::choose_n_farthest_points(kernel, boost::irange(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 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$ + */ + template + Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon) + : Sparse_rips_complex( + boost::irange(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 + 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 + 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 + 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::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; + + if(d * epsilon <= 2 * lj) + alpha = d / 2; + else if(d * epsilon <= li + lj && (epsilon >= 1 || d * epsilon <= lj * (1 + 1 / (1 - epsilon)))) + alpha = d - lj / epsilon; + else continue; + + add_edge(pi, pj, alpha, graph_); + } + } + + Graph graph_; +}; + +} // namespace rips_complex + +} // namespace Gudhi + +#endif // SPARSE_RIPS_COMPLEX_H_ -- cgit v1.2.3 From 342a72ccb89f42109e709ae087d01fa6dcf98e39 Mon Sep 17 00:00:00 2001 From: glisse Date: Fri, 26 Jan 2018 10:22:30 +0000 Subject: doc for epsilon git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/sparserips-glisse@3162 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 02ec0665d288e7e68959aa188cd291924a7d1c8f --- src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) (limited to 'src/Rips_complex/include') diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index c5378b6e..e1e3a951 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -46,7 +46,7 @@ namespace rips_complex { * \ingroup rips_complex * * \details - * This class is used to construct a sparse \f$(1+\epsilon)\f$-approximation of `Rips_complex`. + * This class is used to construct a sparse \f$(1+\epsilon)\f$-approximation of `Rips_complex`, i.e. a filtered simplicial complex that is multiplicatively \f$(1+\epsilon)\f$-interleaved with the Rips filtration. * * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex. */ @@ -65,10 +65,12 @@ 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 (1+epsilon) is the desired approximation factor. epsilon must be positive. * */ template Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) { + GUDHI_CHECK(epsilon > 0, "epsilon must be positive"); std::vector sorted_points; std::vector params; auto dist_fun = [&](Vertex_handle i, Vertex_handle j){return distance(points[i], points[j]);}; @@ -83,6 +85,7 @@ class Sparse_rips_complex { * `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 (1+epsilon) is the desired approximation factor. epsilon must be positive. */ template Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon) -- cgit v1.2.3 From 258439dfbe0d109b1ab26e9b60370cafdc78ac3e Mon Sep 17 00:00:00 2001 From: glisse Date: Tue, 30 Jan 2018 16:31:13 +0000 Subject: Double all filtration values. It was using the Cech convention, not the Rips, so the Rips approximation was not close to the Rips, which is likely to confuse users. git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/sparserips-glisse@3188 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 7b948c956f39dc75931cb00a64516d08ff42e435 --- src/Rips_complex/doc/Intro_rips_complex.h | 4 +++- src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 5 +++-- src/Rips_complex/utilities/README | 12 ++++++------ 3 files changed, 12 insertions(+), 9 deletions(-) (limited to 'src/Rips_complex/include') diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index 54afad66..a9cfc0ab 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -83,7 +83,9 @@ namespace rips_complex { * diagrams are at log-bottleneck distance at most ε). * * The sparse Rips filtration was introduced by Don Sheehy \cite - * sheehy13linear. We are using the version from \cite buchet16efficient. + * sheehy13linear. We are using the version from \cite buchet16efficient + * (except that we multiply all filtration values by 2, to match the usual + * Rips complex). * A more intuitive presentation of the idea is available in \cite * cavanna15geometric, and in a video \cite cavanna15visualizing. * diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index e1e3a951..b8850d00 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -152,10 +152,11 @@ class Sparse_rips_complex { 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 / 2; + alpha = d; else if(d * epsilon <= li + lj && (epsilon >= 1 || d * epsilon <= lj * (1 + 1 / (1 - epsilon)))) - alpha = d - lj / epsilon; + alpha = (d - lj / epsilon) * 2; else continue; add_edge(pi, pj, alpha, graph_); diff --git a/src/Rips_complex/utilities/README b/src/Rips_complex/utilities/README index 40fc58ba..390e31f8 100644 --- a/src/Rips_complex/utilities/README +++ b/src/Rips_complex/utilities/README @@ -97,14 +97,14 @@ where `dim` is the dimension of the homological feature, `birth` and `death` are * `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. **Example with Z/2Z coefficients** -`rips_persistence ../../data/points/tore3D_1307.off -e .5 -m .1 -d 3 -p 2` +`rips_persistence ../../data/points/tore3D_1307.off -e .5 -m .2 -d 3 -p 2` -outputs: +may output (there is randomness involved): ``` -The complex contains 16742034 simplices +The complex contains 17002096 simplices and has dimension 3 2 0 0 inf -2 1 0.0521735 0.645093 -2 2 0.0691677 0.289148 -2 1 0.0491747 0.265232 +2 1 0.104347 1.29343 +2 2 0.138335 0.577743 +2 1 0.0983494 0.529708 ``` -- cgit v1.2.3 From c6957d5b3d8e410b6104a848c793d633fba8990a Mon Sep 17 00:00:00 2001 From: glisse Date: Thu, 1 Feb 2018 22:00:19 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/sparserips-glisse@3208 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 2e0fd8691e04dd1088673910258017802370d265 --- src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/Rips_complex/include') diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index b8850d00..44b73930 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -79,7 +79,7 @@ class Sparse_rips_complex { compute_sparse_graph(sorted_points, params, dist_fun, epsilon); } - /** \brief Rips_complex constructor from a distance matrix. + /** \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 -- cgit v1.2.3 From bae3210f6fb5ef1c498ef669c067b21cd6ddf953 Mon Sep 17 00:00:00 2001 From: glisse Date: Fri, 2 Feb 2018 13:31:23 +0000 Subject: epsilon git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/sparserips-glisse@3209 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: c7cac815867c94821d13d8dd972e6505c43b04ba --- src/Rips_complex/doc/Intro_rips_complex.h | 6 ++++-- src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 8 ++++---- 2 files changed, 8 insertions(+), 6 deletions(-) (limited to 'src/Rips_complex/include') diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index f9663315..909735cd 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -82,12 +82,14 @@ namespace rips_complex { * filtered simplicial complex (linear size, with constants that depend on * ε and the doubling dimension of the space) that is * \f$(1+O(\epsilon))-\f$interleaved with it (in particular, their persistence - * diagrams are at log-bottleneck distance at most ε). + * diagrams are at log-bottleneck distance at most \f$O(\epsilon)\f$). * * The sparse Rips filtration was introduced by Don Sheehy * \cite sheehy13linear. We are using the version from \cite buchet16efficient * (except that we multiply all filtration values by 2, to match the usual - * Rips complex). + * Rips complex), which proves a + * \f$\frac{1+\epsilon}{1-\epsilon}\f$-interleaving, although in practice the + * error is usually smaller. * A more intuitive presentation of the idea is available in * \cite cavanna15geometric, and in a video \cite cavanna15visualizing. * diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 44b73930..503a783a 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -46,7 +46,7 @@ namespace rips_complex { * \ingroup rips_complex * * \details - * This class is used to construct a sparse \f$(1+\epsilon)\f$-approximation of `Rips_complex`, i.e. a filtered simplicial complex that is multiplicatively \f$(1+\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. */ @@ -64,8 +64,8 @@ class Sparse_rips_complex { /** \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 (1+epsilon) is the desired approximation factor. epsilon must be positive. + * @param[in] distance Distance function that returns a `Filtration_value` from 2 given points. + * @param[in] epsilon Approximation parameter. epsilon must be positive. * */ template @@ -85,7 +85,7 @@ class Sparse_rips_complex { * `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 (1+epsilon) is the desired approximation factor. epsilon must be positive. + * @param[in] epsilon Approximation parameter. epsilon must be positive. */ template Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon) -- cgit v1.2.3