diff options
| -rw-r--r-- | biblio/bibliography.bib | 88 | ||||
| -rw-r--r-- | src/Rips_complex/doc/Intro_rips_complex.h | 9 | ||||
| -rw-r--r-- | src/Rips_complex/example/CMakeLists.txt | 6 | ||||
| -rw-r--r-- | src/Rips_complex/example/example_sparse_rips.cpp | 16 | ||||
| -rw-r--r-- | src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 236 | ||||
| -rw-r--r-- | src/Rips_complex/utilities/CMakeLists.txt | 2 | ||||
| -rw-r--r-- | src/Rips_complex/utilities/sparse_rips_persistence.cpp | 58 | ||||
| -rw-r--r-- | src/cython/doc/cubical_complex_user.rst | 2 |
8 files changed, 202 insertions, 215 deletions
diff --git a/biblio/bibliography.bib b/biblio/bibliography.bib index f2ff368d..16f43d6f 100644 --- a/biblio/bibliography.bib +++ b/biblio/bibliography.bib @@ -1,12 +1,12 @@ @inproceedings{gudhilibrary_ICMS14, - author = {Cl\'ement Maria and Jean-Daniel Boissonnat and + author = {Cl\'ement Maria and Jean-Daniel Boissonnat and Marc Glisse and Mariette Yvinec}, title = {The {G}udhi Library: Simplicial Complexes and Persistent Homology}, booktitle = {ICMS}, year = {2014}, } -@article{Carriere17c, +@article{Carriere17c, author = {Carri\`ere, Mathieu and Michel, Bertrand and Oudot, Steve}, title = {{Statistical Analysis and Parameter Selection for Mapper}}, journal = {CoRR}, @@ -20,7 +20,7 @@ booktitle = {Proceedings of the Twenty-ninth Annual Symposium on Computational Geometry}, year = {2013}, pages = {107--116}, -} +} @article{Carriere16, title={{Structure and Stability of the 1-Dimensional Mapper}}, @@ -70,7 +70,7 @@ language={English} publisher = {ACM Press/Addison-Wesley Publishing Co.}, address = {New York, NY, USA}, keywords = {level of detail, mutiresolution modeling, non-manifold, pair contraction, surface simplification}, -} +} @inproceedings{Lindstrom, @@ -85,7 +85,7 @@ language={English} numpages = {8}, publisher = {IEEE Computer Society Press}, address = {Los Alamitos, CA, USA}, -} +} @article{boissonnatmariacamalgorithmica, @@ -201,13 +201,13 @@ language={English}, bibsource = {DBLP, http://dblp.uni-trier.de} } @article{RS62, - author={J. B. Rosser and L. Schoenfeld}, - title={Approximate Formulas for some Functions of Prime Numbers}, - journal= {ijm}, - volume= 6, - year= 1962, - pages={64-94}, - mrnumber={25 \#1139} + author={J. B. Rosser and L. Schoenfeld}, + title={Approximate Formulas for some Functions of Prime Numbers}, + journal= {ijm}, + volume= 6, + year= 1962, + pages={64-94}, + mrnumber={25 \#1139} } @article{DBLP:journals/siamcomp/Furer09, author = {Martin F{\"u}rer}, @@ -410,7 +410,7 @@ url="http://dx.doi.org/10.1007/s00454-013-9557-2" edition = {2}, publisher = {Cambridge University Press}, address = {New York, NY, USA}, -} +} @book{DBLP:books/aw/AhoHU74, author = {Alfred V. Aho and John E. Hopcroft and @@ -508,7 +508,7 @@ proceedings{DBLP:conf/issac/1996, pages = {501--509}, numpages = {9}, keywords = {clique complexes, data structure, edge contraction, flag complexes, high dimensions, homotopy equivalence, shape reconstruction, shape simplification, simplicial complexes, vietoris-rips complexes}, -} +} @@ -571,7 +571,7 @@ note = "http://gmplib.org/", publisher = {ACM}, address = {New York, NY, USA}, keywords = {discrete and computational geometry, embedding theorem, persistent homology}, -} +} @article{Fredman:1984:SST:828.1884, author = {Fredman, Michael L. and Koml\'{o}s, J\'{a}nos and Szemer{\'e}di, Endre}, title = {Storing a Sparse Table with {O}(1) Worst Case Access Time}, @@ -589,7 +589,7 @@ note = "http://gmplib.org/", acmid = {1884}, publisher = {ACM}, address = {New York, NY, USA}, -} +} %------------------------------------------------------------------ @@ -739,7 +739,7 @@ inproceedings{DBLP:conf/soda/1997, publisher = {Kluwer Academic Publishers}, address = {Hingham, MA, USA}, keywords = {Filtration, Klein bottle, Manifold, Natural images, Persistent homology, Topology}, -} +} article{DBLP:journals/ijcv/LeePM03, author = {Ann B. Lee and @@ -767,7 +767,7 @@ pages={234115}, year={2010}, abstract={Understanding energy landscapes is a major challenge in chemistry and biology. Although a wide variety of methods have been invented and applied to this problem, very little is understood about the actual mathematical structures underlying such landscapes. Perhaps the most general assumption is the idea that energy landscapes are low-dimensional manifolds embedded in high-dimensional Euclidean space. While this is a very mild assumption, we have discovered an example of an energy landscape which is nonmanifold, demonstrating previously unknown mathematical complexity. The example occurs in the energy landscape of cyclo-octane, which was found to have the structure of a reducible algebraic variety, composed of the union of a sphere and a Klein bottle, intersecting in two rings.} } - + book{hatcher2002algebraic, title={Algebraic Topology}, @@ -803,7 +803,7 @@ book{hatcher2002algebraic, isbn = {0070131511}, edition = {2nd}, publisher = {McGraw-Hill Higher Education}, -} +} @article{DBLP:journals/focm/CarlssonS10, author = {Gunnar E. Carlsson and @@ -1017,7 +1017,7 @@ language={English} journal = {Journal of Symbolic Computation.}, year = {2016} } - + @ARTICLE{Fasy_Kim_Lecci_Maria_tda, author = {B. Fasy and J. Kim and F. Lecci and C. Maria}, title = {Introduction to the R package TDA.}, @@ -1064,42 +1064,42 @@ language={English} } @ARTICLE{Carriere_Oudot_Ovsjanikov_top_signatures_3d, - author = {M. Carrière and S. Oudot and M. Ovsjanikov}, + author = {M. Carri\`ere and S. Oudot and M. Ovsjanikov}, title = {Stable Topological Signatures for Points on 3D Shapes.}, journal = {Proc. Sympos. on Geometry Processing}, year = {2015} } - + @article{buchet16efficient, - Title = {Efficient and Robust Persistent Homology for Measures}, - Author = {Micka\"{e}l Buchet and Fr\'{e}d\'{e}ric Chazal and Steve Y. Oudot and Donald Sheehy}, - Booktitle = {Computational Geometry: Theory and Applications}, - Volume = {58} - Pages = {70--96} - Year = {2016} + title = {Efficient and Robust Persistent Homology for Measures}, + author = {Micka\"{e}l Buchet and Fr\'{e}d\'{e}ric Chazal and Steve Y. Oudot and Donald Sheehy}, + booktitle = {Computational Geometry: Theory and Applications}, + volume = {58}, + pages = {70--96}, + year = {2016} } @inproceedings{cavanna15geometric, - Author = {Nicholas J. Cavanna and Mahmoodreza Jahanseir and Donald R. Sheehy}, - Booktitle = {Proceedings of the Canadian Conference on Computational Geometry}, - Title = {A Geometric Perspective on Sparse Filtrations}, - Year = {2015} + author = {Nicholas J. Cavanna and Mahmoodreza Jahanseir and Donald R. Sheehy}, + booktitle = {Proceedings of the Canadian Conference on Computational Geometry}, + title = {A Geometric Perspective on Sparse Filtrations}, + year = {2015} } @inproceedings{cavanna15visualizing, - Author = {Nicholas J. Cavanna and Mahmoodreza Jahanseir and Donald R. Sheehy}, - Booktitle = {Proceedings of the 31st International Symposium on Computational Geometry}, - Title = {Visualizing Sparse Filtrations}, - Year = {2015}, - doi = {10.4230/LIPIcs.SOCG.2015.23} + author = {Nicholas J. Cavanna and Mahmoodreza Jahanseir and Donald R. Sheehy}, + booktitle = {Proceedings of the 31st International Symposium on Computational Geometry}, + title = {Visualizing Sparse Filtrations}, + year = {2015}, + doi = {10.4230/LIPIcs.SOCG.2015.23} } @article{sheehy13linear, - Title = {Linear-Size Approximations to the {V}ietoris-{R}ips Filtration}, - Author = {Donald R. Sheehy}, - Journal = {Discrete \& Computational Geometry}, - Volume = {49}, - Number = {4}, - Pages = {778--796}, - Year = {2013} + title = {Linear-Size Approximations to the {V}ietoris-{R}ips Filtration}, + author = {Donald R. Sheehy}, + journal = {Discrete \& Computational Geometry}, + volume = {49}, + number = {4}, + pages = {778--796}, + year = {2013} } diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index da6bd8d9..5a551e60 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -36,7 +36,8 @@ namespace rips_complex { * \section ripsdefinition Rips complex definition * * The Vietoris-Rips complex - * <a target="_blank" href="https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex">(Wikipedia)</a> is an abstract simplicial complex + * <a target="_blank" href="https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex">(Wikipedia)</a> + * is an abstract simplicial complex * defined on a finite metric space, where each simplex corresponds to a subset * of point whose diameter is smaller that some given threshold. * Varying the threshold, we can also see the Rips complex as a filtration of @@ -73,7 +74,8 @@ namespace rips_complex { * * 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"> - * rips_persistence_step_by_step.cpp</a> example, where the constructions of the graph and the Simplex_tree are more detailed. + * rips_persistence_step_by_step.cpp</a> example, where the constructions of the graph and + * the Simplex_tree are more detailed. * * \section sparserips Sparse Rips complex * @@ -147,7 +149,8 @@ namespace rips_complex { * * \subsection sparseripspointscloudexample Example of a sparse Rips from a point cloud * - * This example builds the full sparse Rips of a set of 2D Euclidean points, then prints some minimal information about the 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 * diff --git a/src/Rips_complex/example/CMakeLists.txt b/src/Rips_complex/example/CMakeLists.txt index 90c2d9f7..af86636b 100644 --- a/src/Rips_complex/example/CMakeLists.txt +++ b/src/Rips_complex/example/CMakeLists.txt @@ -11,7 +11,7 @@ add_executable ( Rips_complex_example_one_skeleton_from_distance_matrix example_ add_executable ( Rips_complex_example_from_csv_distance_matrix example_rips_complex_from_csv_distance_matrix_file.cpp ) -# Point cloud +# Sparse rips from points add_executable ( Rips_complex_example_sparse example_sparse_rips.cpp ) # Correlation matrix @@ -22,14 +22,16 @@ if (TBB_FOUND) target_link_libraries(Rips_complex_example_one_skeleton_from_points ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_one_skeleton_from_distance_matrix ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_from_csv_distance_matrix ${TBB_LIBRARIES}) - target_link_libraries(Rips_complex_example_one_skeleton_rips_from_correlation_matrix ${TBB_LIBRARIES}) target_link_libraries(Rips_complex_example_sparse ${TBB_LIBRARIES}) + target_link_libraries(Rips_complex_example_one_skeleton_rips_from_correlation_matrix ${TBB_LIBRARIES}) endif() add_test(NAME Rips_complex_example_one_skeleton_from_points COMMAND $<TARGET_FILE:Rips_complex_example_one_skeleton_from_points>) add_test(NAME Rips_complex_example_one_skeleton_from_distance_matrix COMMAND $<TARGET_FILE:Rips_complex_example_one_skeleton_from_distance_matrix>) +add_test(NAME Rips_complex_example_sparse + COMMAND $<TARGET_FILE:Rips_complex_example_sparse>) add_test(NAME Rips_complex_example_one_skeleton_rips_from_correlation_matrix COMMAND $<TARGET_FILE:Rips_complex_example_one_skeleton_rips_from_correlation_matrix>) diff --git a/src/Rips_complex/example/example_sparse_rips.cpp b/src/Rips_complex/example/example_sparse_rips.cpp index 94e345ea..1c95b48c 100644 --- a/src/Rips_complex/example/example_sparse_rips.cpp +++ b/src/Rips_complex/example/example_sparse_rips.cpp @@ -11,19 +11,12 @@ int main() { using Filtration_value = Simplex_tree::Filtration_value; using Sparse_rips = Gudhi::rips_complex::Sparse_rips_complex<Filtration_value>; - Point points[] = { - {1.0, 1.0}, - {7.0, 0.0}, - {4.0, 6.0}, - {9.0, 6.0}, - {0.0, 14.0}, - {2.0, 19.0}, - {9.0, 17.0}}; + Point points[] = {{1.0, 1.0}, {7.0, 0.0}, {4.0, 6.0}, {9.0, 6.0}, {0.0, 14.0}, {2.0, 19.0}, {9.0, 17.0}}; // ---------------------------------------------------------------------------- // Init from Euclidean points // ---------------------------------------------------------------------------- - double epsilon = 2; // very rough, no guarantees + double epsilon = 2; // very rough, no guarantees Sparse_rips sparse_rips(points, Gudhi::Euclidean_distance(), epsilon); Simplex_tree stree; @@ -32,7 +25,6 @@ int main() { // ---------------------------------------------------------------------------- // Display information about the complex // ---------------------------------------------------------------------------- - std::cout << "Sparse Rips complex is of dimension " << stree.dimension() << - " - " << stree.num_simplices() << " simplices - " << - stree.num_vertices() << " vertices." << std::endl; + std::cout << "Sparse Rips complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() + << " simplices - " << stree.num_vertices() << " vertices." << std::endl; } 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 diff --git a/src/Rips_complex/utilities/CMakeLists.txt b/src/Rips_complex/utilities/CMakeLists.txt index c8d389a0..deb73ff0 100644 --- a/src/Rips_complex/utilities/CMakeLists.txt +++ b/src/Rips_complex/utilities/CMakeLists.txt @@ -27,7 +27,7 @@ add_test(NAME Rips_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:ri add_test(NAME Rips_complex_utility_from_rips_correlation_matrix COMMAND $<TARGET_FILE:rips_correlation_matrix_persistence> "${CMAKE_SOURCE_DIR}/data/correlation_matrix/lower_triangular_correlation_matrix.csv" "-c" "0.3" "-d" "3" "-p" "3" "-m" "0") add_test(NAME Sparse_rips_complex_utility_on_tore_3D COMMAND $<TARGET_FILE:sparse_rips_persistence> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-e" "0.5" "-m" "0.2" "-d" "3" "-p" "2") + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-e" "0.5" "-m" "0.2" "-d" "3" "-p" "2") install(TARGETS rips_distance_matrix_persistence DESTINATION bin) install(TARGETS rips_persistence DESTINATION bin) diff --git a/src/Rips_complex/utilities/sparse_rips_persistence.cpp b/src/Rips_complex/utilities/sparse_rips_persistence.cpp index dd0b2592..d4bae3ba 100644 --- a/src/Rips_complex/utilities/sparse_rips_persistence.cpp +++ b/src/Rips_complex/utilities/sparse_rips_persistence.cpp @@ -1,5 +1,5 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ +/* 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, Clément Maria @@ -36,19 +36,14 @@ using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persis using Filtration_value = Simplex_tree::Filtration_value; using Sparse_rips = Gudhi::rips_complex::Sparse_rips_complex<Filtration_value>; using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; -using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; using Point = std::vector<double>; using Points_off_reader = Gudhi::Points_off_reader<Point>; -void program_options(int argc, char * argv[] - , std::string & off_file_points - , std::string & filediag - , double & epsilon - , int & dim_max - , int & p - , Filtration_value & min_persistence); +void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, double& epsilon, + int& dim_max, int& p, Filtration_value& min_persistence); -int main(int argc, char * argv[]) { +int main(int argc, char* argv[]) { std::string off_file_points; std::string filediag; double epsilon; @@ -90,32 +85,26 @@ int main(int argc, char * argv[]) { return 0; } -void program_options(int argc, char * argv[] - , std::string & off_file_points - , std::string & filediag - , double & epsilon - , int & dim_max - , int & p - , Filtration_value & min_persistence) { +void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, double& epsilon, + int& dim_max, int& p, Filtration_value& min_persistence) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); - hidden.add_options() - ("input-file", po::value<std::string>(&off_file_points), - "Name of an OFF file containing a point set.\n"); + hidden.add_options()("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a point set.\n"); po::options_description visible("Allowed options", 100); - visible.add_options() - ("help,h", "produce help message") - ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::cout") - ("approximation,e", po::value<double>(&epsilon)->default_value(.5), - "Epsilon, where the sparse Rips complex is a (1+epsilon)-approximation of the Rips complex.") - ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), - "Maximal dimension of the Rips complex we want to compute.") - ("field-charac,p", po::value<int>(&p)->default_value(11), - "Characteristic p of the coefficient field Z/pZ for computing homology.") - ("min-persistence,m", po::value<Filtration_value>(&min_persistence), - "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "approximation,e", po::value<double>(&epsilon)->default_value(.5), + "Epsilon, where the sparse Rips complex is a (1+epsilon)-approximation of the Rips complex.")( + "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute.")( + "field-charac,p", po::value<int>(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value<Filtration_value>(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length " + "intervals"); po::positional_options_description pos; pos.add("input-file", 1); @@ -124,8 +113,7 @@ void program_options(int argc, char * argv[] all.add(visible).add(hidden); po::variables_map vm; - po::store(po::command_line_parser(argc, argv). - options(all).positional(pos).run(), vm); + po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm); po::notify(vm); if (vm.count("help") || !vm.count("input-file")) { diff --git a/src/cython/doc/cubical_complex_user.rst b/src/cython/doc/cubical_complex_user.rst index 34598f02..dd82ad93 100644 --- a/src/cython/doc/cubical_complex_user.rst +++ b/src/cython/doc/cubical_complex_user.rst @@ -152,6 +152,6 @@ End user programs are available in cython/example/ folder. Bibliography ============ -.. bibliography:: ../../bibliography.bib +.. bibliography:: ../../biblio/bibliography.bib :filter: docnames :style: unsrt |