diff options
Diffstat (limited to 'src')
52 files changed, 1551 insertions, 228 deletions
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 63c6675c..91305032 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -34,6 +34,7 @@ #include <CGAL/Epick_d.h> #include <CGAL/Spatial_sort_traits_adapter_d.h> #include <CGAL/property_map.h> // for CGAL::Identity_property_map +#include <CGAL/NT_converter.h> #include <iostream> #include <vector> @@ -323,8 +324,9 @@ class Alpha_complex { if (f_simplex_dim > 0) { // squared_radius function initialization Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object(); + CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv; - alpha_complex_filtration = squared_radius(pointVector.begin(), pointVector.end()); + alpha_complex_filtration = cv(squared_radius(pointVector.begin(), pointVector.end())); } complex.assign_filtration(f_simplex, alpha_complex_filtration); #ifdef DEBUG_TRACES diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md index aace85d3..ede749a9 100644 --- a/src/Alpha_complex/utilities/alphacomplex.md +++ b/src/Alpha_complex/utilities/alphacomplex.md @@ -1,3 +1,13 @@ +---
+layout: page
+title: "Alpha complex"
+meta_title: "Alpha complex"
+teaser: ""
+permalink: /alphacomplex/
+---
+{::comment}
+Leave the lines above as it is required by the web site generator 'Jekyll'
+{:/comment}
# Alpha complex #
diff --git a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h index 969daba6..770eb55f 100644 --- a/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h +++ b/src/Bitmap_cubical_complex/include/gudhi/Bitmap_cubical_complex.h @@ -383,7 +383,7 @@ class Bitmap_cubical_complex : public T { std::vector<std::size_t> bdry = this->get_boundary_of_a_cell(sh); if (globalDbg) { std::cerr << "std::pair<Simplex_handle, Simplex_handle> endpoints( Simplex_handle sh )\n"; - std::cerr << "bdry.size() : " << bdry.size() << std::endl; + std::cerr << "bdry.size() : " << bdry.size() << "\n"; } // this method returns two first elements from the boundary of sh. if (bdry.size() < 2) diff --git a/src/Bitmap_cubical_complex/utilities/cubicalcomplex.md b/src/Bitmap_cubical_complex/utilities/cubicalcomplex.md index 6e1b2578..a1bb9007 100644 --- a/src/Bitmap_cubical_complex/utilities/cubicalcomplex.md +++ b/src/Bitmap_cubical_complex/utilities/cubicalcomplex.md @@ -1,3 +1,13 @@ +---
+layout: page
+title: "Cubical complex"
+meta_title: "Cubical complex"
+teaser: ""
+permalink: /cubicalcomplex/
+---
+{::comment}
+Leave the lines above as it is required by the web site generator 'Jekyll'
+{:/comment}
# Cubical complex#
diff --git a/src/Bottleneck_distance/utilities/bottleneckdistance.md b/src/Bottleneck_distance/utilities/bottleneckdistance.md index 526f5822..f2749acc 100644 --- a/src/Bottleneck_distance/utilities/bottleneckdistance.md +++ b/src/Bottleneck_distance/utilities/bottleneckdistance.md @@ -1,3 +1,13 @@ +---
+layout: page
+title: "Bottleneck distance"
+meta_title: "Bottleneck distance"
+teaser: ""
+permalink: /bottleneckdistance/
+---
+{::comment}
+Leave the lines above as it is required by the web site generator 'Jekyll'
+{:/comment}
# Bottleneck distance #
diff --git a/src/GudhUI/model/Model.h b/src/GudhUI/model/Model.h index fc284cc6..072d1185 100644 --- a/src/GudhUI/model/Model.h +++ b/src/GudhUI/model/Model.h @@ -74,7 +74,7 @@ class CGAL_geometric_flag_complex_wrapper { // std::cout << "size:" << vertices.size() << std::endl; for (std::size_t i = 0; i < vertices.size(); ++i) for (std::size_t j = i + 1; j < vertices.size(); ++j) - complex_.add_edge(Vertex_handle(vertices[i]), Vertex_handle(vertices[j])); + complex_.add_edge_without_blockers(Vertex_handle(vertices[i]), Vertex_handle(vertices[j])); } } diff --git a/src/GudhUI/utils/Critical_points.h b/src/GudhUI/utils/Critical_points.h index 2a18e079..e7b9ef31 100644 --- a/src/GudhUI/utils/Critical_points.h +++ b/src/GudhUI/utils/Critical_points.h @@ -79,7 +79,7 @@ template<typename SkBlComplex> class Critical_points { unsigned pos = 0; for (Edge e : edges) { std::cout << "edge " << pos++ << "/" << edges.size() << "\n"; - auto eh = filled_complex_.add_edge(e.first, e.second); + auto eh = filled_complex_.add_edge_without_blockers(e.first, e.second); int is_contractible(is_link_reducible(eh)); switch (is_contractible) { diff --git a/src/GudhUI/utils/K_nearest_builder.h b/src/GudhUI/utils/K_nearest_builder.h index 7be0a4f4..4000a331 100644 --- a/src/GudhUI/utils/K_nearest_builder.h +++ b/src/GudhUI/utils/K_nearest_builder.h @@ -81,7 +81,7 @@ template<typename SkBlComplex> class K_nearest_builder { for (auto it = ++search.begin(); it != search.end(); ++it) { Vertex_handle q(std::get<1>(it->first)); if (p != q && complex_.contains_vertex(p) && complex_.contains_vertex(q)) - complex_.add_edge(p, q); + complex_.add_edge_without_blockers(p, q); } } } diff --git a/src/GudhUI/utils/Rips_builder.h b/src/GudhUI/utils/Rips_builder.h index b22f4db6..59b4bee2 100644 --- a/src/GudhUI/utils/Rips_builder.h +++ b/src/GudhUI/utils/Rips_builder.h @@ -60,7 +60,7 @@ template<typename SkBlComplex> class Rips_builder { std::cout.flush(); for (auto q = p; ++q != vertices.end(); /**/) if (squared_eucl_distance(complex_.point(*p), complex_.point(*q)) < 4 * alpha * alpha) - complex_.add_edge(*p, *q); + complex_.add_edge_without_blockers(*p, *q); } std::cout << std::endl; } diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index 40ff7a4a..f5b67be6 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -23,6 +23,11 @@ #ifndef GIC_H_ #define GIC_H_ +#ifdef GUDHI_USE_TBB +#include <tbb/parallel_for.h> +#include <tbb/mutex.h> +#endif + #include <gudhi/Debug_utils.h> #include <gudhi/graph_simplicial_complex.h> #include <gudhi/reader_utils.h> @@ -99,9 +104,8 @@ class Cover_complex { int data_dimension; // dimension of input data. int n; // number of points. - std::map<int, double> func; // function used to compute the output simplicial complex. - std::map<int, double> - func_color; // function used to compute the colors of the nodes of the output simplicial complex. + std::vector<double> func; // function used to compute the output simplicial complex. + std::vector<double> func_color; // function used to compute the colors of the nodes of the output simplicial complex. bool functional_cover = false; // whether we use a cover with preimages of a function or not. Graph one_skeleton_OFF; // one-skeleton given by the input OFF file (if it exists). @@ -114,8 +118,8 @@ class Cover_complex { Persistence_diagram PD; std::vector<double> distribution; - std::map<int, std::vector<int> > - cover; // function associating to each data point its vectors of cover elements to which it belongs. + std::vector<std::vector<int> > + cover; // function associating to each data point the vector of cover elements to which it belongs. std::map<int, std::vector<int> > cover_back; // inverse of cover, in order to get the data points associated to a specific cover element. std::map<int, double> cover_std; // standard function (induced by func) used to compute the extended persistence @@ -140,8 +144,8 @@ class Cover_complex { // Point comparator struct Less { - Less(std::map<int, double> func) { Fct = func; } - std::map<int, double> Fct; + Less(std::vector<double> func) { Fct = func; } + std::vector<double> Fct; bool operator()(int a, int b) { if (Fct[a] == Fct[b]) return a < b; @@ -276,6 +280,7 @@ class Cover_complex { point_cloud.emplace_back(point.begin(), point.begin() + data_dimension); boost::add_vertex(one_skeleton_OFF); vertices.push_back(boost::add_vertex(one_skeleton)); + std::vector<int> dummy; dummy.clear(); cover.push_back(dummy); i++; } } @@ -430,17 +435,29 @@ class Cover_complex { if (distances.size() == 0) compute_pairwise_distances(distance); - // #pragma omp parallel for - for (int i = 0; i < N; i++) { - SampleWithoutReplacement(n, m, samples); - double hausdorff_dist = 0; - for (int j = 0; j < n; j++) { - double mj = distances[j][samples[0]]; - for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]); - hausdorff_dist = std::max(hausdorff_dist, mj); + #ifdef GUDHI_USE_TBB + tbb::parallel_for(0, N, [&](int i){ + SampleWithoutReplacement(n, m, samples); + double hausdorff_dist = 0; + for (int j = 0; j < n; j++) { + double mj = distances[j][samples[0]]; + for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]); + hausdorff_dist = std::max(hausdorff_dist, mj); + } + delta += hausdorff_dist / N; + }); + #else + for (int i = 0; i < N; i++) { + SampleWithoutReplacement(n, m, samples); + double hausdorff_dist = 0; + for (int j = 0; j < n; j++) { + double mj = distances[j][samples[0]]; + for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]); + hausdorff_dist = std::max(hausdorff_dist, mj); + } + delta += hausdorff_dist / N; } - delta += hausdorff_dist / N; - } + #endif if (verbose) std::cout << "delta = " << delta << std::endl; set_graph_from_rips(delta, distance); @@ -465,7 +482,7 @@ class Cover_complex { while (std::getline(input, line)) { std::stringstream stream(line); stream >> f; - func.emplace(i, f); + func.push_back(f); i++; } functional_cover = true; @@ -479,7 +496,7 @@ class Cover_complex { * */ void set_function_from_coordinate(int k) { - for (int i = 0; i < n; i++) func.emplace(i, point_cloud[i][k]); + for (int i = 0; i < n; i++) func.push_back(point_cloud[i][k]); functional_cover = true; cover_name = "coordinate " + std::to_string(k); } @@ -492,7 +509,7 @@ class Cover_complex { */ template <class InputRange> void set_function_from_range(InputRange const& function) { - for (int i = 0; i < n; i++) func.emplace(i, function[i]); + for (int i = 0; i < n; i++) func.push_back(function[i]); functional_cover = true; } @@ -710,37 +727,69 @@ class Cover_complex { funcstd[i] = 0.5 * (u + v); } - if (verbose) std::cout << "Computing connected components..." << std::endl; - // #pragma omp parallel for - for (int i = 0; i < res; i++) { - // Compute connected components - Graph G = one_skeleton.create_subgraph(); - int num = preimages[i].size(); - std::vector<int> component(num); - for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G); - boost::connected_components(G, &component[0]); - int max = 0; - - // For each point in preimage - for (int j = 0; j < num; j++) { - // Update number of components in preimage - if (component[j] > max) max = component[j]; - - // Identify component with Cantor polynomial N^2 -> N - int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2; - - // Update covers - cover[preimages[i][j]].push_back(identifier); - cover_back[identifier].push_back(preimages[i][j]); - cover_fct[identifier] = i; - cover_std[identifier] = funcstd[i]; - cover_color[identifier].second += func_color[preimages[i][j]]; - cover_color[identifier].first += 1; - } + #ifdef GUDHI_USE_TBB + if (verbose) std::cout << "Computing connected components (parallelized)..." << std::endl; + tbb::parallel_for(0, res, [&](int i){ + // Compute connected components + Graph G = one_skeleton.create_subgraph(); + int num = preimages[i].size(); + std::vector<int> component(num); + for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G); + boost::connected_components(G, &component[0]); + int max = 0; + + // For each point in preimage + for (int j = 0; j < num; j++) { + // Update number of components in preimage + if (component[j] > max) max = component[j]; + + // Identify component with Cantor polynomial N^2 -> N + int identifier = ((i + component[j])*(i + component[j]) + 3 * i + component[j]) / 2; + + // Update covers + cover[preimages[i][j]].push_back(identifier); + cover_back[identifier].push_back(preimages[i][j]); + cover_fct[identifier] = i; + cover_std[identifier] = funcstd[i]; + cover_color[identifier].second += func_color[preimages[i][j]]; + cover_color[identifier].first += 1; + } - // Maximal dimension is total number of connected components - id += max + 1; - } + // Maximal dimension is total number of connected components + id += max + 1; + }); + #else + if (verbose) std::cout << "Computing connected components..." << std::endl; + for (int i = 0; i < res; i++) { + // Compute connected components + Graph G = one_skeleton.create_subgraph(); + int num = preimages[i].size(); + std::vector<int> component(num); + for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G); + boost::connected_components(G, &component[0]); + int max = 0; + + // For each point in preimage + for (int j = 0; j < num; j++) { + // Update number of components in preimage + if (component[j] > max) max = component[j]; + + // Identify component with Cantor polynomial N^2 -> N + int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2; + + // Update covers + cover[preimages[i][j]].push_back(identifier); + cover_back[identifier].push_back(preimages[i][j]); + cover_fct[identifier] = i; + cover_std[identifier] = funcstd[i]; + cover_color[identifier].second += func_color[preimages[i][j]]; + cover_color[identifier].first += 1; + } + + // Maximal dimension is total number of connected components + id += max + 1; + } + #endif maximal_dim = id - 1; for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) @@ -803,24 +852,46 @@ class Cover_complex { for (int j = 0; j < n; j++) mindist[j] = std::numeric_limits<double>::max(); // Compute the geodesic distances to subsamples with Dijkstra - // #pragma omp parallel for - for (int i = 0; i < m; i++) { - if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl; - int seed = voronoi_subsamples[i]; - std::vector<double> dmap(n); - boost::dijkstra_shortest_paths( - one_skeleton, vertices[seed], - boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index))); - - for (int j = 0; j < n; j++) - if (mindist[j] > dmap[j]) { - mindist[j] = dmap[j]; - if (cover[j].size() == 0) - cover[j].push_back(i); - else - cover[j][0] = i; - } - } + #ifdef GUDHI_USE_TBB + if (verbose) std::cout << "Computing geodesic distances (parallelized)..." << std::endl; + tbb::mutex coverMutex; tbb::mutex mindistMutex; + tbb::parallel_for(0, m, [&](int i){ + int seed = voronoi_subsamples[i]; + std::vector<double> dmap(n); + boost::dijkstra_shortest_paths( + one_skeleton, vertices[seed], + boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index))); + + coverMutex.lock(); mindistMutex.lock(); + for (int j = 0; j < n; j++) + if (mindist[j] > dmap[j]) { + mindist[j] = dmap[j]; + if (cover[j].size() == 0) + cover[j].push_back(i); + else + cover[j][0] = i; + } + coverMutex.unlock(); mindistMutex.unlock(); + }); + #else + for (int i = 0; i < m; i++) { + if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl; + int seed = voronoi_subsamples[i]; + std::vector<double> dmap(n); + boost::dijkstra_shortest_paths( + one_skeleton, vertices[seed], + boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index))); + + for (int j = 0; j < n; j++) + if (mindist[j] > dmap[j]) { + mindist[j] = dmap[j]; + if (cover[j].size() == 0) + cover[j].push_back(i); + else + cover[j][0] = i; + } + } + #endif for (int i = 0; i < n; i++) { cover_back[cover[i][0]].push_back(i); @@ -860,7 +931,7 @@ class Cover_complex { while (std::getline(input, line)) { std::stringstream stream(line); stream >> f; - func_color.emplace(i, f); + func_color.push_back(f); i++; } color_name = color_file_name; @@ -873,7 +944,7 @@ class Cover_complex { * */ void set_color_from_coordinate(int k = 0) { - for (int i = 0; i < n; i++) func_color[i] = point_cloud[i][k]; + for (int i = 0; i < n; i++) func_color.push_back(point_cloud[i][k]); color_name = "coordinate "; color_name.append(std::to_string(k)); } @@ -885,7 +956,7 @@ class Cover_complex { * */ void set_color_from_vector(std::vector<double> color) { - for (unsigned int i = 0; i < color.size(); i++) func_color[i] = color[i]; + for (unsigned int i = 0; i < color.size(); i++) func_color.push_back(color[i]); } public: // Create a .dot file that can be compiled with neato to produce a .pdf file. @@ -1039,45 +1110,29 @@ class Cover_complex { minf = std::min(minf, it->second); } + // Build filtration for (auto const& simplex : simplices) { - // Add a simplex and a cone on it - std::vector<int> splx = simplex; - splx.push_back(-2); - st.insert_simplex_and_subfaces(splx); + std::vector<int> splx = simplex; splx.push_back(-2); + st.insert_simplex_and_subfaces(splx, -3); } - // Build filtration - for (auto simplex : st.complex_simplex_range()) { - double filta = std::numeric_limits<double>::lowest(); - double filts = filta; - bool ascending = true; - for (auto vertex : st.simplex_vertex_range(simplex)) { - if (vertex == -2) { - ascending = false; - continue; - } - filta = std::max(-2 + (cover_std[vertex] - minf) / (maxf - minf), filta); - filts = std::max(2 - (cover_std[vertex] - minf) / (maxf - minf), filts); - } - if (ascending) - st.assign_filtration(simplex, filta); - else - st.assign_filtration(simplex, filts); + for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) { + int vertex = it->first; float val = it->second; + int vert[] = {vertex}; int edge[] = {vertex, -2}; + st.assign_filtration(st.find(vert), -2 + (val - minf)/(maxf - minf)); + st.assign_filtration(st.find(edge), 2 - (val - minf)/(maxf - minf)); } - int magic[] = {-2}; - st.assign_filtration(st.find(magic), -3); + st.make_filtration_non_decreasing(); // Compute PD - st.initialize_filtration(); - Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(st); - pcoh.init_coefficients(2); + Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(st); pcoh.init_coefficients(2); pcoh.compute_persistent_cohomology(); // Output PD int max_dim = st.dimension(); for (int i = 0; i < max_dim; i++) { std::vector<std::pair<double, double> > bars = pcoh.intervals_in_dimension(i); - int num_bars = bars.size(); + int num_bars = bars.size(); if(i == 0) num_bars -= 1; if(verbose) std::cout << num_bars << " interval(s) in dimension " << i << ":" << std::endl; for (int j = 0; j < num_bars; j++) { double birth = bars[j].first; @@ -1206,8 +1261,7 @@ class Cover_complex { } if (type == "Nerve") { - for(auto& simplex : cover) - simplices.push_back(simplex.second); + for(int i = 0; i < n; i++) simplices.push_back(cover[i]); std::sort(simplices.begin(), simplices.end()); std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end()); simplices.resize(std::distance(simplices.begin(), it)); diff --git a/src/Nerve_GIC/test/test_GIC.cpp b/src/Nerve_GIC/test/test_GIC.cpp index d633753c..e3067d35 100644 --- a/src/Nerve_GIC/test/test_GIC.cpp +++ b/src/Nerve_GIC/test/test_GIC.cpp @@ -39,6 +39,7 @@ BOOST_AUTO_TEST_CASE(check_nerve) { N.set_type("Nerve"); std::string cloud_file_name("data/cloud"); N.read_point_cloud(cloud_file_name); + N.set_color_from_coordinate(); std::string graph_file_name("data/graph"); N.set_graph_from_file(graph_file_name); std::string cover_file_name("data/cover"); @@ -58,6 +59,7 @@ BOOST_AUTO_TEST_CASE(check_GIC) { GIC.set_type("GIC"); std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name); + GIC.set_color_from_coordinate(); std::string graph_file_name("data/graph"); GIC.set_graph_from_file(graph_file_name); std::string cover_file_name("data/cover"); @@ -77,6 +79,7 @@ BOOST_AUTO_TEST_CASE(check_voronoiGIC) { GIC.set_type("GIC"); std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name); + GIC.set_color_from_coordinate(); std::string graph_file_name("data/graph"); GIC.set_graph_from_file(graph_file_name); GIC.set_cover_from_Voronoi(Gudhi::Euclidean_distance(), 2); diff --git a/src/Nerve_GIC/utilities/covercomplex.md b/src/Nerve_GIC/utilities/covercomplex.md index f33cb2e0..6d16d16f 100644 --- a/src/Nerve_GIC/utilities/covercomplex.md +++ b/src/Nerve_GIC/utilities/covercomplex.md @@ -1,3 +1,13 @@ +--- +layout: page +title: "Cover complex" +meta_title: "Cover complex" +teaser: "" +permalink: /covercomplex/ +--- +{::comment} +Leave the lines above as it is required by the web site generator 'Jekyll' +{:/comment} # Cover complex # diff --git a/src/Persistence_representations/test/CMakeLists.txt b/src/Persistence_representations/test/CMakeLists.txt index 335a71ef..ca1713c3 100644 --- a/src/Persistence_representations/test/CMakeLists.txt +++ b/src/Persistence_representations/test/CMakeLists.txt @@ -35,6 +35,11 @@ target_link_libraries(Read_persistence_from_file_test_unit ${Boost_UNIT_TEST_FRA gudhi_add_coverage_test(Read_persistence_from_file_test_unit) +add_executable ( kernels_unit kernels.cpp ) +target_link_libraries(kernels_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + +gudhi_add_coverage_test(kernels_unit) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) add_executable (Persistence_intervals_with_distances_test_unit persistence_intervals_with_distances_test.cpp ) target_link_libraries(Persistence_intervals_with_distances_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) diff --git a/src/Persistence_representations/test/kernels.cpp b/src/Persistence_representations/test/kernels.cpp new file mode 100644 index 00000000..9db19123 --- /dev/null +++ b/src/Persistence_representations/test/kernels.cpp @@ -0,0 +1,54 @@ +/* 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): Mathieu Carrière + * + * 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 <http://www.gnu.org/licenses/>. + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "kernel" + +#include <boost/test/unit_test.hpp> +#include <cmath> // float comparison +#include <limits> +#include <string> +#include <vector> +#include <algorithm> // std::max +#include <gudhi/common_persistence_representations.h> +#include <gudhi/Weight_functions.h> +#include <gudhi/Persistence_weighted_gaussian.h> +#include <gudhi/Sliced_Wasserstein.h> +#include <gudhi/distance_functions.h> +#include <gudhi/reader_utils.h> + +using SW = Gudhi::Persistence_representations::Sliced_Wasserstein; +using PWG = Gudhi::Persistence_representations::Persistence_weighted_gaussian; + +BOOST_AUTO_TEST_CASE(check_PWG) { + Persistence_diagram v1, v2; v1.emplace_back(0,1); v2.emplace_back(0,2); + PWG pwg1(v1, 1.0, 1000, Gudhi::Persistence_representations::arctan_weight(1,1)); PWG pwgex1(v1, 1.0, -1, Gudhi::Persistence_representations::arctan_weight(1,1)); + PWG pwg2(v2, 1.0, 1000, Gudhi::Persistence_representations::arctan_weight(1,1)); PWG pwgex2(v2, 1.0, -1, Gudhi::Persistence_representations::arctan_weight(1,1)); + BOOST_CHECK(std::abs(pwg1.compute_scalar_product(pwg2) - pwgex1.compute_scalar_product(pwgex2)) <= 1e-1); +} + +BOOST_AUTO_TEST_CASE(check_SW) { + Persistence_diagram v1, v2; v1.emplace_back(0,1); v2.emplace_back(0,2); + SW sw1(v1, 1.0, 100); SW swex1(v1, 1.0, -1); + SW sw2(v2, 1.0, 100); SW swex2(v2, 1.0, -1); + BOOST_CHECK(std::abs(sw1.compute_scalar_product(sw2) - swex1.compute_scalar_product(swex2)) <= 1e-1); +} diff --git a/src/Persistent_cohomology/concept/FilteredComplex.h b/src/Persistent_cohomology/concept/FilteredComplex.h index c19698df..d6b662e9 100644 --- a/src/Persistent_cohomology/concept/FilteredComplex.h +++ b/src/Persistent_cohomology/concept/FilteredComplex.h @@ -31,7 +31,7 @@ struct FilteredComplex typedef unspecified Simplex_handle; /** \brief Key associated to each simplex. * - * Must be a signed integer type. */ + * Must be an integer type. */ typedef unspecified Simplex_key; /** \brief Type for the value of the filtration function. * @@ -67,8 +67,8 @@ struct FilteredComplex Simplex_key key ( Simplex_handle sh ); /** \brief Returns the simplex that has index idx in the filtration. * - * This is never called on null_key(). */ - Simplex_handle simplex ( Simplex_key idx ); + * This is only called on valid indices. */ + Simplex_handle simplex ( size_t idx ); /** \brief Assign a key to a simplex. */ void assign_key(Simplex_handle sh, Simplex_key key); diff --git a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h index 4dbe82c7..3113a22c 100644 --- a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h +++ b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h @@ -162,6 +162,19 @@ persistence diagram with a family of field coefficients. Rips_complex/rips_distance_matrix_persistence.cpp</a> computes the Rips complex of a distance matrix and outputs its persistence diagram. +The file should contain square or lower triangular distance matrix with semicolons as separators. +The code do not check if it is dealing with a distance matrix. It is the user responsibility to provide a valid input. +Please refer to data/distance_matrix/lower_triangular_distance_matrix.csv for an example of a file. + +\li <a href="_rips_complex_2rips_correlation_matrix_persistence_8cpp-example.html"> +Rips_complex/rips_correlation_matrix_persistence.cpp</a> +computes the Rips complex of a correlation matrix and outputs its persistence diagram. + +Note that no check is performed if the matrix given as the input is a correlation matrix. +It is the user responsibility to ensure that this is the case. The input is to be given either as a square or a lower +triangular matrix. +Please refer to data/correlation_matrix/lower_triangular_correlation_matrix.csv for an example of a file. + \li <a href="_alpha_complex_2alpha_complex_3d_persistence_8cpp-example.html"> Alpha_complex/alpha_complex_3d_persistence.cpp</a> computes the persistent homology with \f$\mathbb{Z}/2\mathbb{Z}\f$ coefficients of the alpha complex on points sampling from an OFF file. diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h index e0a147b3..a8c9afa3 100644 --- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h +++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h @@ -285,7 +285,7 @@ class Persistent_cohomology { } } cpx_->assign_key(sigma, cpx_->null_key()); - } else { // If ku == kv, same connected component: create a 1-cocycle class. + } else if (dim_max_ > 1) { // If ku == kv, same connected component: create a 1-cocycle class. create_cocycle(sigma, coeff_field_.multiplicative_identity(), coeff_field_.characteristic()); } } diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index 8c517516..5a551e60 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -29,28 +29,41 @@ namespace rips_complex { /** \defgroup rips_complex Rips complex * - * \author Clément Maria, Pawel Dlotko, Vincent Rouvreau + * \author Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse * * @{ * * \section ripsdefinition Rips complex definition * - * Rips_complex - * <a target="_blank" href="https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex">(Wikipedia)</a> is a - * one skeleton graph that allows to construct a - * <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a> - * from it. - * The input can be a point cloud with a given distance function, or a distance matrix. - * - * The filtration value of each edge is computed from a user-given distance function, or directly from the distance - * matrix. + * 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 + * 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 + * the \f$(n-1)-\f$dimensional simplex, where the filtration value of each + * simplex is the diameter of the corresponding subset of points. + * + * This filtered complex is most often used as an approximation of the + * ÄŒech complex. After rescaling (Rips using the length of the edges and ÄŒech + * the half-length), they share the same 1-skeleton and are multiplicatively + * 2-interleaved or better. While it is slightly bigger, it is also much + * easier to compute. + * + * The number of simplices in the full Rips complex is exponential in the + * number of vertices, it is thus usually restricted, by excluding all the + * simplices with filtration value larger than some threshold, and keeping only + * the dim_max-skeleton. + * + * In order to build this complex, the algorithm first builds the graph. + * The filtration value of each edge is computed from a user-given distance + * function, or directly read from the distance matrix. + * In a second step, this graph is inserted in a simplicial complex, which then + * gets expanded to a flag complex. * - * All edges that have a filtration value strictly greater than a given threshold value are not inserted into - * the complex. + * The input can be given as a range of points and a distance function, or as a + * distance matrix. * - * When creating a simplicial complex from this one skeleton graph, Rips inserts the one skeleton graph into the data - * structure, and then expands the simplicial complex when required. - * * Vertex name correspond to the index of the point in the given range (aka. the point cloud). * * \image html "rips_complex_representation.png" "Rips-complex one skeleton graph representation" @@ -61,7 +74,36 @@ 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 graph construction over the Simplex_tree is 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 + * + * Even truncated in filtration value and dimension, the Rips complex remains + * quite large. However, it is possible to approximate it by a much smaller + * 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 \f$O(\epsilon)\f$). + * + * The sparse Rips filtration was introduced by Don Sheehy + * \cite sheehy13linear. We are using the version described in + * \cite buchet16efficient (except that we multiply all filtration values + * by 2, to match the usual 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. + * + * The interface of `Sparse_rips_complex` is similar to the one for the usual + * `Rips_complex`, except that one has to specify the approximation factor, and + * there is no option to limit the maximum filtration value (the way the + * approximation is done means that larger filtration values are much cheaper + * to handle than low filtration values, so the gain would be too small). + * + * Theoretical guarantees are only available for \f$\epsilon<1\f$. The + * construction accepts larger values of ε, and the size of the complex + * keeps decreasing, but there is no guarantee on the quality of the result. * * \section ripspointsdistance Point cloud and distance function * @@ -104,6 +146,24 @@ namespace rips_complex { * * \include Rips_complex/full_skeleton_rips_for_doc.txt * + * + * \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. + * + * \include Rips_complex/example_sparse_rips.cpp + * + * When launching: + * + * \code $> ./Rips_complex_example_sparse + * \endcode + * + * the program output may be (the exact output varies from one run to the next): + * + * \code Sparse Rips complex is of dimension 2 - 19 simplices - 7 vertices. + * \endcode + * * * * \section ripsdistancematrix Distance matrix @@ -146,6 +206,33 @@ namespace rips_complex { * * \include Rips_complex/full_skeleton_rips_for_doc.txt * + * + * \section ripscorrelationematrix Correlation matrix + * + * Analogously to the case of distance matrix, Rips complexes can be also constructed based on correlation matrix. + * Given a correlation matrix M, comportment-wise 1-M is a distance matrix. + * This example builds the one skeleton graph from the given corelation matrix and threshold value. + * Then it creates a `Simplex_tree` with it. + * + * Then, it is asked to display information about the simplicial complex. + * + * \include Rips_complex/example_one_skeleton_rips_from_correlation_matrix.cpp + * + * When launching: + * + * \code $> ./example_one_skeleton_from_correlation_matrix + * \endcode + * + * the program output is: + * + * \include Rips_complex/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. + * + * @warning As persistence diagrams points will be under the diagonal, bottleneck distance and persistence graphical + * tool will not work properly, this is a known issue. + * */ /** @} */ // end defgroup rips_complex diff --git a/src/Rips_complex/example/CMakeLists.txt b/src/Rips_complex/example/CMakeLists.txt index 2940f164..af86636b 100644 --- a/src/Rips_complex/example/CMakeLists.txt +++ b/src/Rips_complex/example/CMakeLists.txt @@ -11,17 +11,29 @@ 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 ) +# Sparse rips from points +add_executable ( Rips_complex_example_sparse example_sparse_rips.cpp ) + +# Correlation matrix +add_executable ( Rips_complex_example_one_skeleton_rips_from_correlation_matrix example_one_skeleton_rips_from_correlation_matrix.cpp ) + if (TBB_FOUND) target_link_libraries(Rips_complex_example_from_off ${TBB_LIBRARIES}) 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_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>) add_test(NAME Rips_complex_example_from_off_doc_12_1 COMMAND $<TARGET_FILE:Rips_complex_example_from_off> "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "12.0" "1" "${CMAKE_CURRENT_BINARY_DIR}/ripsoffreader_result_12_1.txt") @@ -57,3 +69,4 @@ install(TARGETS Rips_complex_example_from_off DESTINATION bin) install(TARGETS Rips_complex_example_one_skeleton_from_points DESTINATION bin) install(TARGETS Rips_complex_example_one_skeleton_from_distance_matrix DESTINATION bin) install(TARGETS Rips_complex_example_from_csv_distance_matrix DESTINATION bin) +install(TARGETS Rips_complex_example_one_skeleton_rips_from_correlation_matrix DESTINATION bin) diff --git a/src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp b/src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp new file mode 100644 index 00000000..1343f24d --- /dev/null +++ b/src/Rips_complex/example/example_one_skeleton_rips_from_correlation_matrix.cpp @@ -0,0 +1,81 @@ +#include <gudhi/Rips_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> + +#include <iostream> +#include <string> +#include <vector> +#include <limits> // for std::numeric_limits + +int main() { + // Type definitions + using Simplex_tree = Gudhi::Simplex_tree<>; + using Filtration_value = Simplex_tree::Filtration_value; + using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; + using Distance_matrix = std::vector<std::vector<Filtration_value>>; + + // User defined correlation matrix is: + // |1 0.06 0.23 0.01 0.89| + // |0.06 1 0.74 0.01 0.61| + // |0.23 0.74 1 0.72 0.03| + // |0.01 0.01 0.72 1 0.7 | + // |0.89 0.61 0.03 0.7 1 | + + Distance_matrix correlations; + correlations.push_back({}); + correlations.push_back({0.06}); + correlations.push_back({0.23, 0.74}); + correlations.push_back({0.01, 0.01, 0.72}); + correlations.push_back({0.89, 0.61, 0.03, 0.7}); + + // ---------------------------------------------------------------------------- + // Convert correlation matrix to a distance matrix: + // ---------------------------------------------------------------------------- + double threshold = 0; + for (size_t i = 0; i != correlations.size(); ++i) { + for (size_t j = 0; j != correlations[i].size(); ++j) { + // Here we check if our data comes from corelation matrix. + if ((correlations[i][j] < -1) || (correlations[i][j] > 1)) { + std::cerr << "The input matrix is not a correlation matrix. The program will now terminate.\n"; + throw "The input matrix is not a correlation matrix. The program will now terminate.\n"; + } + correlations[i][j] = 1 - correlations[i][j]; + // Here we make sure that we will get the treshold value equal to maximal + // distance in the matrix. + if (correlations[i][j] > threshold) threshold = correlations[i][j]; + } + } + + //----------------------------------------------------------------------------- + // Now the correlation matrix is a distance matrix and can be processed further. + //----------------------------------------------------------------------------- + Distance_matrix distances = correlations; + + Rips_complex rips_complex_from_points(distances, threshold); + + Simplex_tree stree; + rips_complex_from_points.create_complex(stree, 1); + // ---------------------------------------------------------------------------- + // Display information about the one skeleton Rips complex. Note that + // the filtration displayed here comes from the distance matrix computed + // above, which is 1 - initial correlation matrix. Only this way, we obtain + // a complex with filtration. If a correlation matrix is used instead, we would + // have a reverse filtration (i.e. filtration of boundary of each simplex S + // is greater or equal to the filtration of S). + // ---------------------------------------------------------------------------- + std::cout << "Rips complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() << " simplices - " + << stree.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on Rips complex simplices in the filtration order, with [filtration value]:" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + std::cout << " ( "; + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " + << "[" << stree.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + + return 0; +} diff --git a/src/Rips_complex/example/example_sparse_rips.cpp b/src/Rips_complex/example/example_sparse_rips.cpp new file mode 100644 index 00000000..1c95b48c --- /dev/null +++ b/src/Rips_complex/example/example_sparse_rips.cpp @@ -0,0 +1,30 @@ +#include <gudhi/Sparse_rips_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> + +#include <iostream> +#include <vector> + +int main() { + using Point = std::vector<double>; + using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; + 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}}; + + // ---------------------------------------------------------------------------- + // Init from Euclidean points + // ---------------------------------------------------------------------------- + double epsilon = 2; // very rough, no guarantees + Sparse_rips sparse_rips(points, Gudhi::Euclidean_distance(), epsilon); + + Simplex_tree stree; + sparse_rips.create_complex(stree, 10); + + // ---------------------------------------------------------------------------- + // 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; +} diff --git a/src/Rips_complex/example/one_skeleton_rips_from_correlation_matrix_for_doc.txt b/src/Rips_complex/example/one_skeleton_rips_from_correlation_matrix_for_doc.txt new file mode 100644 index 00000000..640d7083 --- /dev/null +++ b/src/Rips_complex/example/one_skeleton_rips_from_correlation_matrix_for_doc.txt @@ -0,0 +1,17 @@ +Rips complex is of dimension 1 - 15 simplices - 5 vertices. +Iterator on Rips complex simplices in the filtration order, with [filtration value]: + ( 0 ) -> [0] + ( 1 ) -> [0] + ( 2 ) -> [0] + ( 3 ) -> [0] + ( 4 ) -> [0] + ( 4 0 ) -> [0.11] + ( 2 1 ) -> [0.26] + ( 3 2 ) -> [0.28] + ( 4 3 ) -> [0.3] + ( 4 1 ) -> [0.39] + ( 2 0 ) -> [0.77] + ( 1 0 ) -> [0.94] + ( 4 2 ) -> [0.97] + ( 3 0 ) -> [0.99] + ( 3 1 ) -> [0.99] 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..1a9d6ebb --- /dev/null +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -0,0 +1,175 @@ +/* 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef SPARSE_RIPS_COMPLEX_H_ +#define SPARSE_RIPS_COMPLEX_H_ + +#include <gudhi/Debug_utils.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/choose_n_farthest_points.h> + +#include <boost/graph/adjacency_list.hpp> +#include <boost/range/metafunctions.hpp> + +#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. + +/** + * \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. + * + * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex. + */ +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_); + } + } + + Graph graph_; +}; + +} // namespace rips_complex + +} // namespace Gudhi + +#endif // SPARSE_RIPS_COMPLEX_H_ diff --git a/src/Rips_complex/test/test_rips_complex.cpp b/src/Rips_complex/test/test_rips_complex.cpp index 89afbc25..4e7b79d2 100644 --- a/src/Rips_complex/test/test_rips_complex.cpp +++ b/src/Rips_complex/test/test_rips_complex.cpp @@ -31,6 +31,7 @@ #include <algorithm> // std::max #include <gudhi/Rips_complex.h> +#include <gudhi/Sparse_rips_complex.h> // to construct Rips_complex from a OFF file of points #include <gudhi/Points_off_io.h> #include <gudhi/Simplex_tree.h> @@ -43,6 +44,7 @@ using Point = std::vector<double>; using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; using Rips_complex = Gudhi::rips_complex::Rips_complex<Simplex_tree::Filtration_value>; +using Sparse_rips_complex = Gudhi::rips_complex::Sparse_rips_complex<Simplex_tree::Filtration_value>; using Distance_matrix = std::vector<std::vector<Filtration_value>>; BOOST_AUTO_TEST_CASE(RIPS_DOC_OFF_file) { @@ -230,6 +232,71 @@ BOOST_AUTO_TEST_CASE(Rips_complex_from_points) { } } +BOOST_AUTO_TEST_CASE(Sparse_rips_complex_from_points) { + // This is a clone of the test above + // ---------------------------------------------------------------------------- + // Init of a list of points + // ---------------------------------------------------------------------------- + Vector_of_points points; + std::vector<double> coords = { 0.0, 0.0, 0.0, 1.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 0.0, 0.0, 1.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 0.0, 1.0, 0.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 1.0, 0.0, 0.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + + // ---------------------------------------------------------------------------- + // Init of a Rips complex from the list of points + // ---------------------------------------------------------------------------- + // .001 is small enough that we get a deterministic result matching the exact Rips + Sparse_rips_complex sparse_rips(points, Custom_square_euclidean_distance(), .001); + + std::cout << "========== Sparse_rips_complex_from_points ==========" << std::endl; + Simplex_tree st; + const int DIMENSION = 3; + sparse_rips.create_complex(st, DIMENSION); + + // Another way to check num_simplices + std::cout << "Iterator on Rips complex simplices in the filtration order, with [filtration value]:" << std::endl; + int num_simplices = 0; + for (auto f_simplex : st.filtration_simplex_range()) { + num_simplices++; + std::cout << " ( "; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << st.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + BOOST_CHECK(num_simplices == 15); + std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl; + BOOST_CHECK(st.num_simplices() == 15); + + std::cout << "st.dimension()=" << st.dimension() << std::endl; + BOOST_CHECK(st.dimension() == DIMENSION); + std::cout << "st.num_vertices()=" << st.num_vertices() << std::endl; + BOOST_CHECK(st.num_vertices() == 4); + + for (auto f_simplex : st.filtration_simplex_range()) { + std::cout << "dimension(" << st.dimension(f_simplex) << ") - f = " << st.filtration(f_simplex) << std::endl; + switch (st.dimension(f_simplex)) { + case 0: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.0); + break; + case 1: + case 2: + case 3: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 2.0); + break; + default: + BOOST_CHECK(false); // Shall not happen + break; + } + } +} + BOOST_AUTO_TEST_CASE(Rips_doc_csv_file) { // ---------------------------------------------------------------------------- // diff --git a/src/Rips_complex/utilities/CMakeLists.txt b/src/Rips_complex/utilities/CMakeLists.txt index baa571fa..deb73ff0 100644 --- a/src/Rips_complex/utilities/CMakeLists.txt +++ b/src/Rips_complex/utilities/CMakeLists.txt @@ -7,15 +7,29 @@ target_link_libraries(rips_distance_matrix_persistence ${Boost_PROGRAM_OPTIONS_L add_executable(rips_persistence rips_persistence.cpp) target_link_libraries(rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY}) +add_executable(rips_correlation_matrix_persistence rips_correlation_matrix_persistence.cpp) +target_link_libraries(rips_correlation_matrix_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +add_executable(sparse_rips_persistence sparse_rips_persistence.cpp) +target_link_libraries(sparse_rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY}) + if (TBB_FOUND) target_link_libraries(rips_distance_matrix_persistence ${TBB_LIBRARIES}) target_link_libraries(rips_persistence ${TBB_LIBRARIES}) + target_link_libraries(rips_correlation_matrix_persistence ${TBB_LIBRARIES}) + target_link_libraries(sparse_rips_persistence ${TBB_LIBRARIES}) endif() add_test(NAME Rips_complex_utility_from_rips_distance_matrix COMMAND $<TARGET_FILE:rips_distance_matrix_persistence> "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" "-r" "1.0" "-d" "3" "-p" "3" "-m" "0") add_test(NAME Rips_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:rips_persistence> "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") +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_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) +install(TARGETS rips_correlation_matrix_persistence DESTINATION bin) +install(TARGETS sparse_rips_persistence DESTINATION bin) diff --git a/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp new file mode 100644 index 00000000..c2082fae --- /dev/null +++ b/src/Rips_complex/utilities/rips_correlation_matrix_persistence.cpp @@ -0,0 +1,170 @@ +/* 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): Pawel Dlotko, Vincent Rouvreau + * + * Copyright (C) 2016 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 <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Rips_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/reader_utils.h> +#include <gudhi/writing_persistence_to_file.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::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 Correlation_matrix = std::vector<std::vector<Filtration_value>>; +using intervals_common = Gudhi::Persistence_interval_common<double, int>; + +void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag, + Filtration_value& correlation_min, int& dim_max, int& p, Filtration_value& min_persistence); + +int main(int argc, char* argv[]) { + std::string csv_matrix_file; + std::string filediag; + Filtration_value correlation_min; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, csv_matrix_file, filediag, correlation_min, dim_max, p, min_persistence); + + Correlation_matrix correlations = + Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_matrix_file); + + Filtration_value threshold = 0; + + // Given a correlation matrix M, we compute component-wise M'[i,j] = 1-M[i,j] to get a distance matrix: + for (size_t i = 0; i != correlations.size(); ++i) { + for (size_t j = 0; j != correlations[i].size(); ++j) { + correlations[i][j] = 1 - correlations[i][j]; + // Here we make sure that the values of corelations lie between -1 and 1. + // If not, we throw an exception. + if ((correlations[i][j] < -1) || (correlations[i][j] > 1)) { + std::cerr << "The input matrix is not a correlation matrix. The program will now terminate. \n"; + throw "The input matrix is not a correlation matrix. The program will now terminate. \n"; + } + if (correlations[i][j] > threshold) threshold = correlations[i][j]; + } + } + + Rips_complex rips_complex_from_file(correlations, threshold); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree simplex_tree; + + rips_complex_from_file.create_complex(simplex_tree, dim_max); + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + // compute persistence + pcoh.compute_persistent_cohomology(min_persistence); + + // invert the persistence diagram. The reason for this procedure is the following: + // The input to the program is a corelation matrix M. When processing it, it is + // turned into 1-M and the obtained persistence intervals are in '1-M' units. + // Below we reverse every (birth,death) pair into (1-birth, 1-death) pair + // so that the input and the output to the program is expressed in the same + // units. + auto pairs = pcoh.get_persistent_pairs(); + std::vector<intervals_common> processed_persistence_intervals; + processed_persistence_intervals.reserve(pairs.size()); + for (auto pair : pairs) { + double birth = 1 - simplex_tree.filtration(get<0>(pair)); + double death = 1 - simplex_tree.filtration(get<1>(pair)); + unsigned dimension = (unsigned)simplex_tree.dimension(get<0>(pair)); + int field = get<2>(pair); + processed_persistence_intervals.push_back(intervals_common(birth, death, dimension, field)); + } + + // sort the processed intervals: + std::sort(processed_persistence_intervals.begin(), processed_persistence_intervals.end()); + + // and write them to a file + if (filediag.empty()) { + write_persistence_intervals_to_stream(processed_persistence_intervals); + } else { + std::ofstream out(filediag); + write_persistence_intervals_to_stream(processed_persistence_intervals, out); + } + return 0; +} + +void program_options(int argc, char* argv[], std::string& csv_matrix_file, std::string& filediag, + Filtration_value& correlation_min, 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>(&csv_matrix_file), + "Name of file containing a corelation matrix. Can be square or lower triangular matrix. Separator is ';'."); + 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")( + "min-edge-corelation,c", po::value<Filtration_value>(&correlation_min)->default_value(0), + "Minimal corelation of an edge for the Rips complex construction.")( + "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); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map 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")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a corelation matrix.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Rips_complex/utilities/ripscomplex.md b/src/Rips_complex/utilities/ripscomplex.md index 4291fae7..84f40cc1 100644 --- a/src/Rips_complex/utilities/ripscomplex.md +++ b/src/Rips_complex/utilities/ripscomplex.md @@ -1,3 +1,13 @@ +--- +layout: page +title: "Rips complex" +meta_title: "Rips complex" +teaser: "" +permalink: /ripscomplex/ +--- +{::comment} +Leave the lines above as it is required by the web site generator 'Jekyll' +{:/comment} # Rips complex # @@ -39,11 +49,63 @@ Same as `rips_persistence` but taking a distance matrix as input. **Usage** -`rips_persistence [options] <CSV input file>` +`rips_distance_matrix_persistence [options] <CSV input file>` where `<CSV input file>` is the path to the file containing a distance matrix. Can be square or lower triangular matrix. Separator is ';'. +The code do not check if it is dealing with a distance matrix. It is the user responsibility to provide a valid input. +Please refer to data/distance_matrix/lower_triangular_distance_matrix.csv for an example of a file. **Example** `rips_distance_matrix_persistence data/distance_matrix/full_square_distance_matrix.csv -r 15 -d 3 -p 3 -m 0` + + +## rips_correlation_matrix_persistence ## + +Same as `rips_distance_matrix_persistence` but taking a correlation matrix as input. + +**Usage** + +`rips_correlation_matrix_persistence [options] <CSV input file>` + +where +`<CSV input file>` is the path to the file containing a correlation matrix. Can be square or lower triangular matrix. Separator is ';'. +Note that no check is performed if the matrix given as the input is a correlation matrix. +It is the user responsibility to ensure that this is the case. +Please refer to data/correlation_matrix/lower_triangular_correlation_matrix.csv for an example of a file. + +**Example** + +`rips_correlation_matrix_persistence data/distance_matrix/full_square_distance_matrix.csv -r 15 -d 3 -p 3 -m 0` + +**Warning** + +As persistence diagrams points will be under the diagonal, bottleneck distance and persistence graphical tool will not work +properly, this is a known issue. + + +## sparse_rips_persistence ## +This program computes the persistent homology with coefficient field *Z/pZ* +of a sparse (1+epsilon)-approximation of the Rips complex defined on a set of input Euclidean points. The output diagram contains one bar per line, written with the convention: + +`p dim birth death` + +where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number). + +**Usage** + +`sparse_rips_persistence [options] <OFF input file>` + +**Allowed options** + +* `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. +* `-e [ --approximation ]` (default = .5) Epsilon, where the sparse Rips complex is a (1+epsilon)-approximation of the Rips complex. +* `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Rips complex we want to compute. +* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-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** + +`sparse_rips_persistence ../../data/points/tore3D_1307.off -e .5 -m .2 -d 3 -p 2` diff --git a/src/Rips_complex/utilities/sparse_rips_persistence.cpp b/src/Rips_complex/utilities/sparse_rips_persistence.cpp new file mode 100644 index 00000000..d4bae3ba --- /dev/null +++ b/src/Rips_complex/utilities/sparse_rips_persistence.cpp @@ -0,0 +1,133 @@ +/* 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 + * + * 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 <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Sparse_rips_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +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 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); + +int main(int argc, char* argv[]) { + std::string off_file_points; + std::string filediag; + double epsilon; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, epsilon, dim_max, p, min_persistence); + + Points_off_reader off_reader(off_file_points); + Sparse_rips sparse_rips(off_reader.get_point_cloud(), Gudhi::Euclidean_distance(), epsilon); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree simplex_tree; + + sparse_rips.create_complex(simplex_tree, dim_max); + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } + + 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) { + 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"); + + 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"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map 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")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a sparse (1+epsilon)-approximation of the Rips complex \ndefined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Witness_complex/utilities/witnesscomplex.md b/src/Witness_complex/utilities/witnesscomplex.md index 2341759b..3be9bc55 100644 --- a/src/Witness_complex/utilities/witnesscomplex.md +++ b/src/Witness_complex/utilities/witnesscomplex.md @@ -1,3 +1,13 @@ +---
+layout: page
+title: "Witness complex"
+meta_title: "Witness complex"
+teaser: ""
+permalink: /witnesscomplex/
+---
+{::comment}
+Leave the lines above as it is required by the web site generator 'Jekyll'
+{:/comment}
# Witness complex #
diff --git a/src/common/doc/header.html b/src/common/doc/header.html index d69b28fa..2f54e68d 100644 --- a/src/common/doc/header.html +++ b/src/common/doc/header.html @@ -9,7 +9,7 @@ <!--BEGIN PROJECT_NAME--><title>$projectname: $title</title><!--END PROJECT_NAME--> <!--BEGIN !PROJECT_NAME--><title>$title</title><!--END !PROJECT_NAME--> <!-- GUDHI website css for header BEGIN --> -<link rel="stylesheet" type="text/css" href="http://pages.saclay.inria.fr/vincent.rouvreau/gudhi/gudhi-doc-2.0.0/assets/css/styles_feeling_responsive.css" /> +<link rel="stylesheet" type="text/css" href="http://gudhi.gforge.inria.fr/assets/css/styles_feeling_responsive.css" /> <!-- GUDHI website css for header END --> <link href="$relpath^tabs.css" rel="stylesheet" type="text/css"/> <script type="text/javascript" src="$relpath^jquery.js"></script> diff --git a/src/common/doc/offline_header.html b/src/common/doc/offline_header.html new file mode 100644 index 00000000..6a02a895 --- /dev/null +++ b/src/common/doc/offline_header.html @@ -0,0 +1,41 @@ +<!-- HTML header for doxygen 1.8.6--> +<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd"> +<!-- GUDHI website : class="no-js" lang="en" is necessary --> +<html xmlns="http://www.w3.org/1999/xhtml" class="no-js" lang="en"> +<head> +<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/> +<meta http-equiv="X-UA-Compatible" content="IE=9"/> +<meta name="generator" content="Doxygen $doxygenversion"/> +<!--BEGIN PROJECT_NAME--><title>$projectname: $title</title><!--END PROJECT_NAME--> +<!--BEGIN !PROJECT_NAME--><title>$title</title><!--END !PROJECT_NAME--> +<!-- GUDHI website css for header END --> +<link href="$relpath^tabs.css" rel="stylesheet" type="text/css"/> +<script type="text/javascript" src="$relpath^jquery.js"></script> +<script type="text/javascript" src="$relpath^dynsections.js"></script> +$treeview +$search +$mathjax +<link href="$relpath^$stylesheet" rel="stylesheet" type="text/css" /> +$extrastylesheet +</head> +<body> + + +<div id="top"><!-- do not remove this div, it is closed by doxygen! --> + +<!--BEGIN TITLEAREA--> +<div id="titlearea"> +<table cellspacing="0" cellpadding="0"> + <tbody> + <tr style="height: 30px;"> + <!--BEGIN DISABLE_INDEX--> + <!--BEGIN SEARCHENGINE--> + <td>$searchbox</td> + <!--END SEARCHENGINE--> + <!--END DISABLE_INDEX--> + </tr> + </tbody> +</table> +</div> +<!--END TITLEAREA--> +<!-- end header part --> diff --git a/src/common/example/CMakeLists.txt b/src/common/example/CMakeLists.txt index afe865d4..1273c699 100644 --- a/src/common/example/CMakeLists.txt +++ b/src/common/example/CMakeLists.txt @@ -3,11 +3,20 @@ project(Common_examples) add_executable ( vector_double_off_reader example_vector_double_points_off_reader.cpp ) target_link_libraries(vector_double_off_reader ${CGAL_LIBRARY}) +file(COPY "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) add_test(NAME Common_example_vector_double_off_reader COMMAND $<TARGET_FILE:vector_double_off_reader> - "${CMAKE_SOURCE_DIR}/data/points/SO3_10000.off") + "alphacomplexdoc.off") install(TARGETS vector_double_off_reader DESTINATION bin) +if (DIFF_PATH) + # Do not forget to copy test results files in current binary dir + file(COPY "vectordoubleoffreader_result.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + + add_test(Common_example_vector_double_off_reader_diff_files ${DIFF_PATH} + ${CMAKE_CURRENT_BINARY_DIR}/vectordoubleoffreader_result.txt ${CMAKE_CURRENT_BINARY_DIR}/alphacomplexdoc.off.txt) +endif() + if(CGAL_FOUND) add_executable ( cgal_3D_off_reader example_CGAL_3D_points_off_reader.cpp ) target_link_libraries(cgal_3D_off_reader ${CGAL_LIBRARY}) diff --git a/src/common/example/example_CGAL_3D_points_off_reader.cpp b/src/common/example/example_CGAL_3D_points_off_reader.cpp index 665b7a29..4658d8d5 100644 --- a/src/common/example/example_CGAL_3D_points_off_reader.cpp +++ b/src/common/example/example_CGAL_3D_points_off_reader.cpp @@ -20,12 +20,12 @@ int main(int argc, char **argv) { usage(argv[0]); } - std::string offInputFile(argv[1]); + std::string off_input_file(argv[1]); // Read the OFF file (input file name given as parameter) and triangulate points - Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile); + Gudhi::Points_3D_off_reader<Point_3> off_reader(off_input_file); // Check the read operation was correct if (!off_reader.is_valid()) { - std::cerr << "Unable to read file " << offInputFile << std::endl; + std::cerr << "Unable to read file " << off_input_file << std::endl; usage(argv[0]); } diff --git a/src/common/example/example_CGAL_points_off_reader.cpp b/src/common/example/example_CGAL_points_off_reader.cpp index 8c6a6b54..f45683a5 100644 --- a/src/common/example/example_CGAL_points_off_reader.cpp +++ b/src/common/example/example_CGAL_points_off_reader.cpp @@ -22,12 +22,12 @@ int main(int argc, char **argv) { usage(argv[0]); } - std::string offInputFile(argv[1]); + std::string off_input_file(argv[1]); // Read the OFF file (input file name given as parameter) and triangulate points - Gudhi::Points_off_reader<Point_d> off_reader(offInputFile); + Gudhi::Points_off_reader<Point_d> off_reader(off_input_file); // Check the read operation was correct if (!off_reader.is_valid()) { - std::cerr << "Unable to read file " << offInputFile << std::endl; + std::cerr << "Unable to read file " << off_input_file << std::endl; usage(argv[0]); } diff --git a/src/common/example/example_vector_double_points_off_reader.cpp b/src/common/example/example_vector_double_points_off_reader.cpp index 8aecb26e..5093da85 100644 --- a/src/common/example/example_vector_double_points_off_reader.cpp +++ b/src/common/example/example_vector_double_points_off_reader.cpp @@ -17,25 +17,27 @@ int main(int argc, char **argv) { usage(argv[0]); } - std::string offInputFile(argv[1]); + std::string off_input_file(argv[1]); // Read the OFF file (input file name given as parameter) and triangulate points - Gudhi::Points_off_reader<Point_d> off_reader(offInputFile); + Gudhi::Points_off_reader<Point_d> off_reader(off_input_file); // Check the read operation was correct if (!off_reader.is_valid()) { - std::cerr << "Unable to read file " << offInputFile << std::endl; + std::cerr << "Unable to read file " << off_input_file << std::endl; usage(argv[0]); } // Retrieve the triangulation std::vector<Point_d> point_cloud = off_reader.get_point_cloud(); + std::ofstream output_file(off_input_file + ".txt"); int n {0}; for (auto point : point_cloud) { - std::cout << "Point[" << n << "] = "; + output_file << "Point[" << n << "] = "; for (std::size_t i {0}; i < point.size(); i++) - std::cout << point[i] << " "; - std::cout << "\n"; + output_file << point[i] << " "; + output_file << "\n"; ++n; } + output_file.close(); return 0; } diff --git a/src/common/example/cgaloffreader_result.txt b/src/common/example/vectordoubleoffreader_result.txt index 1deb8dbd..1deb8dbd 100644 --- a/src/common/example/cgaloffreader_result.txt +++ b/src/common/example/vectordoubleoffreader_result.txt diff --git a/src/common/include/gudhi/Off_reader.h b/src/common/include/gudhi/Off_reader.h index 4fcd2af2..32320e4d 100644 --- a/src/common/include/gudhi/Off_reader.h +++ b/src/common/include/gudhi/Off_reader.h @@ -105,25 +105,26 @@ class Off_reader { bool is_off_file = (line.find("OFF") != std::string::npos); bool is_noff_file = (line.find("nOFF") != std::string::npos); + + if (!is_off_file && !is_noff_file) { std::cerr << line << std::endl; std::cerr << "missing off header\n"; return false; } + if (is_noff_file) { + // Should be on a separate line, but we accept it on the same line as the number of vertices + stream_ >> off_info_.dim; + } else { + off_info_.dim = 3; + } + if (!goto_next_uncomment_line(line)) return false; std::istringstream iss(line); - if ((is_off_file) && (!is_noff_file)) { - off_info_.dim = 3; - if (!(iss >> off_info_.num_vertices >> off_info_.num_faces >> off_info_.num_edges)) { - std::cerr << "incorrect number of vertices/faces/edges\n"; - return false; - } - } else { - if (!(iss >> off_info_.dim >> off_info_.num_vertices >> off_info_.num_faces >> off_info_.num_edges)) { + if (!(iss >> off_info_.num_vertices >> off_info_.num_faces >> off_info_.num_edges)) { std::cerr << "incorrect number of vertices/faces/edges\n"; return false; - } } off_visitor.init(off_info_.dim, off_info_.num_vertices, off_info_.num_faces, off_info_.num_edges); @@ -131,10 +132,12 @@ class Off_reader { } bool goto_next_uncomment_line(std::string& uncomment_line) { - uncomment_line.clear(); - do - std::getline(stream_, uncomment_line); while (uncomment_line[0] == '%'); - return (uncomment_line.size() > 0 && uncomment_line[0] != '%'); + do { + // skip whitespace, including empty lines + if (!std::ifstream::sentry(stream_)) return false; + std::getline(stream_, uncomment_line); + } while (uncomment_line[0] == '#'); + return (bool)stream_; } template<typename OffVisitor> diff --git a/src/common/include/gudhi/Points_off_io.h b/src/common/include/gudhi/Points_off_io.h index 29af8a8a..08f324c6 100644 --- a/src/common/include/gudhi/Points_off_io.h +++ b/src/common/include/gudhi/Points_off_io.h @@ -126,9 +126,9 @@ class Points_off_visitor_reader { * \code $> ./vector_double_off_reader ../../data/points/alphacomplexdoc.off * \endcode * - * the program output is: + * the program outputs a file ../../data/points/alphacomplexdoc.off.txt: * - * \include common/cgaloffreader_result.txt + * \include common/vectordoubleoffreader_result.txt */ template<typename Point_d> class Points_off_reader { diff --git a/src/common/include/gudhi/writing_persistence_to_file.h b/src/common/include/gudhi/writing_persistence_to_file.h new file mode 100644 index 00000000..4c5ce918 --- /dev/null +++ b/src/common/include/gudhi/writing_persistence_to_file.h @@ -0,0 +1,117 @@ +/* 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): Pawel Dlotko + * + * Copyright (C) 2017 Swansea University, UK + * + * 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 <http://www.gnu.org/licenses/>. + */ + +#ifndef WRITING_PERSISTENCE_TO_FILE_H +#define WRITING_PERSISTENCE_TO_FILE_H + +#include <iostream> +#include <string> +#include <limits> + +namespace Gudhi { + +/** +* This is a class to store persistence intervals. Its main purpose is to +* exchange data in between different packages and provide unified way +* of writing a collection of persistence intervals to file. +**/ +template <typename Filtration_type, typename Coefficient_field> +class Persistence_interval_common { + public: + /** + * Constructor taking as an input birth and death of the pair. + **/ + Persistence_interval_common(Filtration_type birth, Filtration_type death) + : birth_(birth), + death_(death), + dimension_(std::numeric_limits<unsigned>::max()), + arith_element_(std::numeric_limits<Coefficient_field>::max()) {} + + /** + * Constructor taking as an input birth, death and dimension of the pair. + **/ + Persistence_interval_common(Filtration_type birth, Filtration_type death, unsigned dim) + : birth_(birth), death_(death), dimension_(dim), arith_element_(std::numeric_limits<Coefficient_field>::max()) {} + + /** +* Constructor taking as an input birth, death, dimension of the pair as well +* as the number p such that this interval is present over Z_p field. +**/ + Persistence_interval_common(Filtration_type birth, Filtration_type death, unsigned dim, Coefficient_field field) + : birth_(birth), death_(death), dimension_(dim), arith_element_(field) {} + + /** + * Operator to compare two persistence pairs. During the comparision all the + * fields: birth, death, dimensiona and arith_element_ are taken into account + * and they all have to be equal for two pairs to be equal. + **/ + bool operator==(const Persistence_interval_common& i2) const { + return ((this->birth_ == i2.birth_) && (this->death_ == i2.death_) && (this->dimension_ == i2.dimension_) && + (this->arith_element_ == i2.arith_element_)); + } + + /** + * Check if two persistence paris are not equal. + **/ + bool operator!=(const Persistence_interval_common& i2) const { return (!((*this) == i2)); } + + /** + * Operator to compare objects of a type Persistence_interval_common. + * One intervals is smaller than the other if it has lower persistence. + * Note that this operator do not take Arith_element into account when doing comparisions. + **/ + bool operator<(const Persistence_interval_common& i2) const { + return fabs(this->death_ - this->birth_) < fabs(i2.death_ - i2.birth_); + } + + friend std::ostream& operator<<(std::ostream& out, const Persistence_interval_common& it) { + if (it.arith_element_ != std::numeric_limits<Coefficient_field>::max()) { + out << it.arith_element_ << " "; + } + if (it.dimension_ != std::numeric_limits<unsigned>::max()) { + out << it.dimension_ << " "; + } + out << it.birth_ << " " << it.death_ << " "; + return out; + } + + private: + Filtration_type birth_; + Filtration_type death_; + unsigned dimension_; + Coefficient_field arith_element_; +}; + +/** + * This function write a vector<Persistence_interval_common> to a stream +**/ +template <typename Persistence_interval_range> +void write_persistence_intervals_to_stream(const Persistence_interval_range& intervals, + std::ostream& out = std::cout) { + for (auto interval : intervals) { + out << interval << "\n"; + } +} + +} + +#endif // WRITING_PERSISTENCE_TO_FILE_H diff --git a/src/common/utilities/pointsetgenerator.md b/src/common/utilities/pointsetgenerator.md index 284715d4..3b23e668 100644 --- a/src/common/utilities/pointsetgenerator.md +++ b/src/common/utilities/pointsetgenerator.md @@ -1,6 +1,16 @@ - - -# common # +--- +layout: page +title: "OFF point set generator" +meta_title: "OFF point set generator" +teaser: "" +permalink: /pointsetgenerator/ +--- +{::comment} +Leave the lines above as it is required by the web site generator 'Jekyll' +{:/comment} + + +# Miscellaneous # ## off_file_from_shape_generator ## diff --git a/src/cython/cython/off_reader.pyx b/src/cython/cython/off_reader.pyx index b6e107ef..266dae2c 100644 --- a/src/cython/cython/off_reader.pyx +++ b/src/cython/cython/off_reader.pyx @@ -46,4 +46,5 @@ def read_off(off_file=''): return read_points_from_OFF_file(str.encode(off_file)) else: print("file " + off_file + " not found.") + return [] diff --git a/src/cython/cython/rips_complex.pyx b/src/cython/cython/rips_complex.pyx index ad9b0a4d..73b154b8 100644 --- a/src/cython/cython/rips_complex.pyx +++ b/src/cython/cython/rips_complex.pyx @@ -34,8 +34,6 @@ __license__ = "GPL v3" cdef extern from "Rips_complex_interface.h" namespace "Gudhi": cdef cppclass Rips_complex_interface "Gudhi::rips_complex::Rips_complex_interface": Rips_complex_interface(vector[vector[double]] values, double threshold, bool euclidean) - # bool from_file is a workaround for cython to find the correct signature - Rips_complex_interface(string file_name, double threshold, bool euclidean, bool from_file) void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, int dim_max) # RipsComplex python interface @@ -49,7 +47,7 @@ cdef class RipsComplex: cdef Rips_complex_interface * thisptr # Fake constructor that does nothing but documenting the constructor - def __init__(self, points=None, off_file='', distance_matrix=None, csv_file='', max_edge_length=float('inf')): + def __init__(self, points=None, distance_matrix=None, max_edge_length=float('inf')): """RipsComplex constructor. :param max_edge_length: Rips value. @@ -60,41 +58,14 @@ cdef class RipsComplex: Or - :param off_file: An OFF file style name. - :type off_file: string - - Or - :param distance_matrix: A distance matrix (full square or lower triangular). :type points: list of list of double - - Or - - :param csv_file: A csv file style name containing a full square or a - lower triangular distance matrix. - :type csv_file: string """ # The real cython constructor - def __cinit__(self, points=None, off_file='', distance_matrix=None, csv_file='', max_edge_length=float('inf')): - if off_file is not '': - if os.path.isfile(off_file): - self.thisptr = new Rips_complex_interface(str.encode(off_file), - max_edge_length, - True, - True) - else: - print("file " + off_file + " not found.") - elif csv_file is not '': - if os.path.isfile(csv_file): - self.thisptr = new Rips_complex_interface(str.encode(csv_file), - max_edge_length, - False, - True) - else: - print("file " + csv_file + " not found.") - elif distance_matrix is not None: + def __cinit__(self, points=None, distance_matrix=None, max_edge_length=float('inf')): + if distance_matrix is not None: self.thisptr = new Rips_complex_interface(distance_matrix, max_edge_length, False) else: if points is None: 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 diff --git a/src/cython/doc/persistence_graphical_tools_user.rst b/src/cython/doc/persistence_graphical_tools_user.rst index 9033331f..a5523d23 100644 --- a/src/cython/doc/persistence_graphical_tools_user.rst +++ b/src/cython/doc/persistence_graphical_tools_user.rst @@ -58,8 +58,8 @@ This function can display the persistence result as a diagram: import gudhi - rips_complex = gudhi.RipsComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/tore3D_1307.off', max_edge_length=0.2) + point_cloud = gudhi.read_off(off_file=gudhi.__root_source_dir__ + '/data/points/tore3D_1307.off') + rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=0.2) simplex_tree = rips_complex.create_simplex_tree(max_dimension=3) diag = simplex_tree.persistence() plt = gudhi.plot_persistence_diagram(diag, band_boot=0.13) @@ -69,8 +69,8 @@ This function can display the persistence result as a diagram: import gudhi - rips_complex = gudhi.RipsComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/tore3D_1307.off', max_edge_length=0.2) + point_cloud = gudhi.read_off(off_file=gudhi.__root_source_dir__ + '/data/points/tore3D_1307.off') + rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=0.2) simplex_tree = rips_complex.create_simplex_tree(max_dimension=3) diag = simplex_tree.persistence() plt = gudhi.plot_persistence_diagram(diag, band_boot=0.13) diff --git a/src/cython/doc/pyplots/diagram_persistence.py b/src/cython/doc/pyplots/diagram_persistence.py index c2fbf801..ac20bf47 100755 --- a/src/cython/doc/pyplots/diagram_persistence.py +++ b/src/cython/doc/pyplots/diagram_persistence.py @@ -1,7 +1,8 @@ import gudhi -rips_complex = gudhi.RipsComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/tore3D_1307.off', max_edge_length=0.2) +point_cloud = gudhi.read_off(off_file=gudhi.__root_source_dir__ + \ + '/data/points/tore3D_1307.off') +rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=0.2) simplex_tree = rips_complex.create_simplex_tree(max_dimension=3) diag = simplex_tree.persistence() plt = gudhi.plot_persistence_diagram(diag, band_boot=0.13) diff --git a/src/cython/doc/rips_complex_user.rst b/src/cython/doc/rips_complex_user.rst index 96ba9944..7738aef0 100644 --- a/src/cython/doc/rips_complex_user.rst +++ b/src/cython/doc/rips_complex_user.rst @@ -101,8 +101,8 @@ Finally, it is asked to display information about the Rips complex. .. testcode:: import gudhi - rips_complex = gudhi.RipsComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/alphacomplexdoc.off', max_edge_length=12.0) + point_cloud = gudhi.read_off(off_file=gudhi.__root_source_dir__ + '/data/points/alphacomplexdoc.off') + rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=12.0) simplex_tree = rips_complex.create_simplex_tree(max_dimension=1) result_str = 'Rips complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ repr(simplex_tree.num_simplices()) + ' simplices - ' + \ @@ -197,7 +197,7 @@ Example from csv file ^^^^^^^^^^^^^^^^^^^^^ This example builds the :doc:`Rips_complex <rips_complex_ref>` from the given -points in an OFF file, and max_edge_length value. +distance matrix in a csv file, and max_edge_length value. Then it creates a :doc:`Simplex_tree <simplex_tree_ref>` with it. Finally, it is asked to display information about the Rips complex. @@ -206,8 +206,9 @@ Finally, it is asked to display information about the Rips complex. .. testcode:: import gudhi - rips_complex = gudhi.RipsComplex(csv_file=gudhi.__root_source_dir__ + \ - '/data/distance_matrix/full_square_distance_matrix.csv', max_edge_length=12.0) + distance_matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file=gudhi.__root_source_dir__ + \ + '/data/distance_matrix/full_square_distance_matrix.csv') + rips_complex = gudhi.RipsComplex(distance_matrix=distance_matrix, max_edge_length=12.0) simplex_tree = rips_complex.create_simplex_tree(max_dimension=1) result_str = 'Rips complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ repr(simplex_tree.num_simplices()) + ' simplices - ' + \ @@ -240,3 +241,72 @@ the program output is: [0, 3] -> 9.43 [4, 6] -> 9.49 [3, 6] -> 11.00 + +Correlation matrix +--------------- + +Example from a correlation matrix +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Analogously to the case of distance matrix, Rips complexes can be also constructed based on correlation matrix. +Given a correlation matrix M, comportment-wise 1-M is a distance matrix. +This example builds the one skeleton graph from the given corelation matrix and threshold value. +Then it creates a :doc:`Simplex_tree <simplex_tree_ref>` with it. + +Finally, it is asked to display information about the simplicial complex. + +.. testcode:: + + import gudhi + import numpy as np + + # User defined correlation matrix is: + # |1 0.06 0.23 0.01 0.89| + # |0.06 1 0.74 0.01 0.61| + # |0.23 0.74 1 0.72 0.03| + # |0.01 0.01 0.72 1 0.7 | + # |0.89 0.61 0.03 0.7 1 | + correlation_matrix=np.array([[1., 0.06, 0.23, 0.01, 0.89], + [0.06, 1., 0.74, 0.01, 0.61], + [0.23, 0.74, 1., 0.72, 0.03], + [0.01, 0.01, 0.72, 1., 0.7], + [0.89, 0.61, 0.03, 0.7, 1.]], float) + + distance_matrix = np.ones((correlation_matrix.shape),float) - correlation_matrix + rips_complex = gudhi.RipsComplex(distance_matrix=distance_matrix, max_edge_length=1.0) + + simplex_tree = rips_complex.create_simplex_tree(max_dimension=1) + result_str = 'Rips complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ + repr(simplex_tree.num_simplices()) + ' simplices - ' + \ + repr(simplex_tree.num_vertices()) + ' vertices.' + print(result_str) + fmt = '%s -> %.2f' + for filtered_value in simplex_tree.get_filtration(): + print(fmt % tuple(filtered_value)) + +When launching (Rips maximal distance between 2 points is 12.0, is expanded +until dimension 1 - one skeleton graph in other words), the output is: + +.. testoutput:: + + Rips complex is of dimension 1 - 15 simplices - 5 vertices. + [0] -> 0.00 + [1] -> 0.00 + [2] -> 0.00 + [3] -> 0.00 + [4] -> 0.00 + [0, 4] -> 0.11 + [1, 2] -> 0.26 + [2, 3] -> 0.28 + [3, 4] -> 0.30 + [1, 4] -> 0.39 + [0, 2] -> 0.77 + [0, 1] -> 0.94 + [2, 4] -> 0.97 + [0, 3] -> 0.99 + [1, 3] -> 0.99 + +.. note:: + As persistence diagrams points will be under the diagonal, + bottleneck distance and persistence graphical tool will not work properly, + this is a known issue. diff --git a/src/cython/example/alpha_rips_persistence_bottleneck_distance.py b/src/cython/example/alpha_rips_persistence_bottleneck_distance.py index ab5fc1e9..386f8457 100755 --- a/src/cython/example/alpha_rips_persistence_bottleneck_distance.py +++ b/src/cython/example/alpha_rips_persistence_bottleneck_distance.py @@ -45,13 +45,14 @@ args = parser.parse_args() with open(args.file, 'r') as f: first_line = f.readline() if (first_line == 'OFF\n') or (first_line == 'nOFF\n'): + point_cloud = gudhi.read_off(off_file=args.file) print("#####################################################################") print("RipsComplex creation from points read in a OFF file") message = "RipsComplex with max_edge_length=" + repr(args.threshold) print(message) - rips_complex = gudhi.RipsComplex(off_file=args.file, + rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=args.threshold) rips_stree = rips_complex.create_simplex_tree(max_dimension=args.max_dimension) @@ -67,7 +68,7 @@ with open(args.file, 'r') as f: message = "AlphaComplex with max_edge_length=" + repr(args.threshold) print(message) - alpha_complex = gudhi.AlphaComplex(off_file=args.file) + alpha_complex = gudhi.AlphaComplex(points=point_cloud) alpha_stree = alpha_complex.create_simplex_tree(max_alpha_square=(args.threshold * args.threshold)) message = "Number of simplices=" + repr(alpha_stree.num_simplices()) diff --git a/src/cython/example/bottleneck_basic_example.py b/src/cython/example/bottleneck_basic_example.py index 31cecb29..a7fa01c1 100755 --- a/src/cython/example/bottleneck_basic_example.py +++ b/src/cython/example/bottleneck_basic_example.py @@ -28,8 +28,6 @@ __author__ = "Francois Godi, Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 INRIA" __license__ = "GPL v3" -import gudhi - diag1 = [[2.7, 3.7],[9.6, 14.],[34.2, 34.974], [3.,float('Inf')]] diag2 = [[2.8, 4.45],[9.5, 14.1],[3.2,float('Inf')]] diff --git a/src/cython/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py b/src/cython/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py new file mode 100755 index 00000000..aa82ef71 --- /dev/null +++ b/src/cython/example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python + +import gudhi +import sys +import argparse + +"""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): Vincent Rouvreau + + Copyright (C) 2017 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 <http://www.gnu.org/licenses/>. +""" + +__author__ = "Vincent Rouvreau" +__copyright__ = "Copyright (C) 2017 INRIA" +__license__ = "GPL v3" + +parser = argparse.ArgumentParser(description='RipsComplex creation from ' + 'a correlation matrix read in a csv file.', + epilog='Example: ' + 'example/rips_complex_diagram_persistence_from_correlation_matrix_file_example.py ' + '-f ../data/correlation_matrix/lower_triangular_correlation_matrix.csv -e 12.0 -d 3' + '- Constructs a Rips complex with the ' + 'correlation matrix from the given csv file.') +parser.add_argument("-f", "--file", type=str, required=True) +parser.add_argument("-c", "--min_edge_correlation", type=float, default=0.5) +parser.add_argument("-d", "--max_dimension", type=int, default=1) +parser.add_argument("-b", "--band_boot", type=float, default=0.) +parser.add_argument('--no-diagram', default=False, action='store_true' , help='Flag for not to display the diagrams') + +args = parser.parse_args() + +if not (-1. < args.min_edge_correlation < 1.): + print("Wrong value of the treshold corelation (should be between -1 and 1).") + sys.exit(1) + +print("#####################################################################") +print("Caution: as persistence diagrams points will be under the diagonal,") +print("bottleneck distance and persistence graphical tool will not work") +print("properly, this is a known issue.") + +print("#####################################################################") +print("RipsComplex creation from correlation matrix read in a csv file") + +message = "RipsComplex with min_edge_correlation=" + repr(args.min_edge_correlation) +print(message) + +correlation_matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file=args.file) +# Given a correlation matrix M, we compute component-wise M'[i,j] = 1-M[i,j] to get a distance matrix: +distance_matrix = [[1.-correlation_matrix[i][j] for j in range(len(correlation_matrix[i]))] for i in range(len(correlation_matrix))] + +rips_complex = gudhi.RipsComplex(distance_matrix=distance_matrix, + max_edge_length=1.-args.min_edge_correlation) +simplex_tree = rips_complex.create_simplex_tree(max_dimension=args.max_dimension) + +message = "Number of simplices=" + repr(simplex_tree.num_simplices()) +print(message) + +diag = simplex_tree.persistence() + +print("betti_numbers()=") +print(simplex_tree.betti_numbers()) + +# invert the persistence diagram +invert_diag = [(diag[pers][0],(1.-diag[pers][1][0], 1.-diag[pers][1][1])) for pers in range(len(diag))] + +if args.no_diagram == False: + pplot = gudhi.plot_persistence_diagram(invert_diag, band_boot=args.band_boot) + pplot.show() diff --git a/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py b/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py index 3baebd17..c8aac240 100755 --- a/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py +++ b/src/cython/example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py @@ -30,12 +30,12 @@ __copyright__ = "Copyright (C) 2016 INRIA" __license__ = "GPL v3" parser = argparse.ArgumentParser(description='RipsComplex creation from ' - 'a distance matrix read in a OFF file.', + 'a distance matrix read in a csv file.', epilog='Example: ' 'example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py ' '-f ../data/distance_matrix/lower_triangular_distance_matrix.csv -e 12.0 -d 3' '- Constructs a Rips complex with the ' - 'points from the given OFF file.') + 'distance matrix from the given csv file.') parser.add_argument("-f", "--file", type=str, required=True) parser.add_argument("-e", "--max_edge_length", type=float, default=0.5) parser.add_argument("-d", "--max_dimension", type=int, default=1) @@ -50,7 +50,8 @@ print("RipsComplex creation from distance matrix read in a csv file") message = "RipsComplex with max_edge_length=" + repr(args.max_edge_length) print(message) -rips_complex = gudhi.RipsComplex(csv_file=args.file, max_edge_length=args.max_edge_length) +distance_matrix = gudhi.read_lower_triangular_matrix_from_csv_file(csv_file=args.file) +rips_complex = gudhi.RipsComplex(distance_matrix=distance_matrix, max_edge_length=args.max_edge_length) simplex_tree = rips_complex.create_simplex_tree(max_dimension=args.max_dimension) message = "Number of simplices=" + repr(simplex_tree.num_simplices()) diff --git a/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py b/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py index 5951eedf..544b68c9 100755 --- a/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py +++ b/src/cython/example/rips_complex_diagram_persistence_from_off_file_example.py @@ -53,7 +53,8 @@ with open(args.file, 'r') as f: message = "RipsComplex with max_edge_length=" + repr(args.max_edge_length) print(message) - rips_complex = gudhi.RipsComplex(off_file=args.file, max_edge_length=args.max_edge_length) + point_cloud = gudhi.read_off(off_file=args.file) + rips_complex = gudhi.RipsComplex(points=point_cloud, max_edge_length=args.max_edge_length) simplex_tree = rips_complex.create_simplex_tree(max_dimension=args.max_dimension) message = "Number of simplices=" + repr(simplex_tree.num_simplices()) diff --git a/src/cython/include/Rips_complex_interface.h b/src/cython/include/Rips_complex_interface.h index 02985727..f26befbc 100644 --- a/src/cython/include/Rips_complex_interface.h +++ b/src/cython/include/Rips_complex_interface.h @@ -25,9 +25,7 @@ #include <gudhi/Simplex_tree.h> #include <gudhi/Rips_complex.h> -#include <gudhi/Points_off_io.h> #include <gudhi/distance_functions.h> -#include <gudhi/reader_utils.h> #include "Simplex_tree_interface.h" @@ -56,21 +54,6 @@ class Rips_complex_interface { } } - Rips_complex_interface(const std::string& file_name, double threshold, bool euclidean, bool from_file = true) { - if (euclidean) { - // Rips construction where file_name is an OFF file - Gudhi::Points_off_reader<Point_d> off_reader(file_name); - rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(off_reader.get_point_cloud(), - threshold, - Gudhi::Euclidean_distance()); - } else { - // Rips construction where values is a distance matrix - Distance_matrix distances = - Gudhi::read_lower_triangular_matrix_from_csv_file<Simplex_tree_interface<>::Filtration_value>(file_name); - rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(distances, threshold); - } - } - ~Rips_complex_interface() { delete rips_complex_; } |