diff options
Diffstat (limited to 'src')
| -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 |
7 files changed, 158 insertions, 171 deletions
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 |