From 1ce5d0d19e13a14e8a67442aec7bc40eae68dc8e Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Sun, 12 Apr 2020 10:33:38 +0200 Subject: Modify interface and utils --- .../include/gudhi/Flag_complex_sparse_matrix.h | 86 ++++++++------ src/Collapse/test/collapse_unit_test.cpp | 129 +++++++++++++-------- ...tance_matrix_edge_collapse_rips_persistence.cpp | 28 ++--- .../point_cloud_edge_collapse_rips_persistence.cpp | 15 +-- 4 files changed, 148 insertions(+), 110 deletions(-) (limited to 'src') diff --git a/src/Collapse/include/gudhi/Flag_complex_sparse_matrix.h b/src/Collapse/include/gudhi/Flag_complex_sparse_matrix.h index d472bf15..e90d284d 100644 --- a/src/Collapse/include/gudhi/Flag_complex_sparse_matrix.h +++ b/src/Collapse/include/gudhi/Flag_complex_sparse_matrix.h @@ -16,6 +16,7 @@ #include #include +#include #include @@ -32,6 +33,7 @@ #include #include + namespace Gudhi { namespace collapse { @@ -47,26 +49,28 @@ namespace collapse { * Algorithm. * * \tparam Vertex_handle type must be a signed integer type. It admits a total order <. - * \tparam Filtration_value type for the value of the filtration function. Must be comparable with <. + * \tparam Filtration type for the value of the filtration function. Must be comparable with <. */ -template +template class Flag_complex_sparse_matrix { private: using Edge = std::pair; // This is an ordered pair, An edge is stored with convention of the first // element being the smaller i.e {2,3} not {3,2}. However this is at the level // of row indices on actual vertex lables - using Filtered_edge = std::pair; + using Filtered_edge = std::pair; using Row_index = std::size_t; using Map_vertex_to_index = std::unordered_map; - using Sparse_row_matrix = Eigen::SparseMatrix; + using Sparse_row_matrix = Eigen::SparseMatrix; using Row_indices_vector = std::vector; public: - using Filtered_sorted_edge_list = std::vector>; + using Filtered_sorted_edge_list = std::vector>; + using Filtration_value = Filtration; + using Proximity_graph = Gudhi::Proximity_graph; private: std::unordered_map row_to_vertex_; @@ -88,9 +92,9 @@ class Flag_complex_sparse_matrix { // Boolean vector to indicate if the index is critical or not. std::vector dominated_edge_indicator_; // domination indicator - //! Stores the Map between vertices_row_to_vertex_ and row indices row_to_vertex_ -> row-index. + //! Stores the Map between verticesrow_to_vertex_ and row indices row_to_vertex_ -> row-index. /*! - So, if the original simplex tree had vertices_ 0,1,4,5
+ So, if the original simplex tree had vertices 0,1,4,5
row_to_vertex_ would store :
\verbatim Values = | 0 | 1 | 4 | 5 | @@ -106,10 +110,10 @@ class Flag_complex_sparse_matrix { */ std::unordered_map vertex_to_row_; - //! Stores the Sparse matrix of Filtration_value values representing the Original Simplicial Complex. + //! Stores the Sparse matrix of Filtration values representing the Original Simplicial Complex. /*! \code - Sparse_row_matrix = Eigen::SparseMatrix ; + Sparse_row_matrix = Eigen::SparseMatrix ; \endcode ; */ @@ -120,7 +124,7 @@ class Flag_complex_sparse_matrix { //! Stores true for dominated rows and false for undominated rows. /*! Initialised to a vector of length equal to the value of the variable rows with all false values. - Subsequent removal of dominated vertices_ is reflected by concerned entries changing to true in this vector. + Subsequent removal of dominated vertices is reflected by concerned entries changing to true in this vector. */ std::vector domination_indicator_; //(domination indicator) @@ -128,9 +132,9 @@ class Flag_complex_sparse_matrix { std::vector f_edge_vector_; - //! Stores the number of vertices_ in the original Simplicial Complex. + //! Stores the number of vertices in the original Simplicial Complex. /*! - This stores the count of vertices_ (which is also the number of rows in the Matrix). + This stores the count of vertices (which is also the number of rows in the Matrix). */ Row_index rows; @@ -202,7 +206,7 @@ class Flag_complex_sparse_matrix { } template - void set_edge_critical(Row_index indx, Filtration_value filt, FilteredEdgeInsertion filtered_edge_insert) { + void set_edge_critical(Row_index indx, Filtration filt, FilteredEdgeInsertion filtered_edge_insert) { #ifdef DEBUG_TRACES std::cout << "The curent index with filtration value " << indx << ", " << filt << " is primary critical" << std::endl; @@ -290,7 +294,7 @@ class Flag_complex_sparse_matrix { return common; } - void insert_vertex(Vertex_handle vertex, Filtration_value filt_val) { + void insert_vertex(Vertex_handle vertex, Filtration filt_val) { auto rw = vertex_to_row_.find(vertex); if (rw == vertex_to_row_.end()) { // Initializing the diagonal element of the adjency matrix corresponding to rw_b. @@ -302,7 +306,7 @@ class Flag_complex_sparse_matrix { } } - void insert_new_edges(Vertex_handle u, Vertex_handle v, Filtration_value filt_val) + void insert_new_edges(Vertex_handle u, Vertex_handle v, Filtration filt_val) { // The edge must not be added before, it should be a new edge. insert_vertex(u, filt_val); @@ -336,19 +340,27 @@ class Flag_complex_sparse_matrix { domination_indicator_ are initialised by init() function which is called at the begining of this.
*/ - Flag_complex_sparse_matrix(const Filtered_sorted_edge_list& edge_t) + template + Flag_complex_sparse_matrix(const DistanceMatrix& distance_matrix) : rows(0) { - for (size_t bgn_idx = 0; bgn_idx < edge_t.size(); bgn_idx++) { - Vertex_handle u = std::get<1>(edge_t[bgn_idx]); - Vertex_handle v = std::get<2>(edge_t[bgn_idx]); - f_edge_vector_.push_back({{u, v}, std::get<0>(edge_t[bgn_idx])}); - vertices_.emplace(u); - vertices_.emplace(v); + Vertex_handle num_vertices = std::distance(std::begin(distance_matrix), std::end(distance_matrix)); + + // This one is not part of the loop + vertices_.emplace(0); + // Only browse the lower part of the distance matrix + for (Vertex_handle line_index = 1; line_index < num_vertices; line_index++) { + for (Vertex_handle row_index = 0; row_index < line_index; row_index++) { +#ifdef DEBUG_TRACES + std::cout << "Insert edge: fn[" << row_index << ", " << line_index << "] = " + << distance_matrix[line_index][row_index] << std::endl; +#endif // DEBUG_TRACES + f_edge_vector_.push_back({{row_index, line_index}, distance_matrix[line_index][row_index]}); + } + vertices_.emplace(line_index); } } - template - Flag_complex_sparse_matrix(const OneSkeletonGraph& one_skeleton_graph) + Flag_complex_sparse_matrix(const Proximity_graph& one_skeleton_graph) : rows(0) { // Insert all vertices_ for (auto v_it = boost::vertices(one_skeleton_graph); v_it.first != v_it.second; ++v_it.first) { @@ -362,6 +374,18 @@ class Flag_complex_sparse_matrix { Vertex_handle v = target(edge, one_skeleton_graph); f_edge_vector_.push_back({{u, v}, boost::get(Gudhi::edge_filtration_t(), one_skeleton_graph, edge)}); } + } + + // Performs edge collapse in a decreasing sequence of the filtration value. + template + void filtered_edge_collapse(FilteredEdgeInsertion filtered_edge_insert) { + Row_index endIdx = 0; + + u_set_removed_redges_.clear(); + u_set_dominated_redges_.clear(); + critical_edge_indicator_.clear(); + + std::cout << "Sort it - " << f_edge_vector_.size() << std::endl; // Sort edges auto sort_by_filtration = [](const Filtered_edge& edge_a, const Filtered_edge& edge_b) -> bool { @@ -373,17 +397,9 @@ class Flag_complex_sparse_matrix { #else std::stable_sort(f_edge_vector_.begin(), f_edge_vector_.end(), sort_by_filtration); #endif - } - - // Performs edge collapse in a decreasing sequence of the filtration value. - template - void filtered_edge_collapse(FilteredEdgeInsertion filtered_edge_insert) { - Row_index endIdx = 0; - - u_set_removed_redges_.clear(); - u_set_dominated_redges_.clear(); - critical_edge_indicator_.clear(); + std::cout << "Sorted" << std::endl; + std::cout << vertices_.size() << std::endl; // Initializing sparse_row_adjacency_matrix_, This is a row-major sparse matrix. sparse_row_adjacency_matrix_ = Sparse_row_matrix(vertices_.size(), vertices_.size()); @@ -392,7 +408,7 @@ class Flag_complex_sparse_matrix { Edge edge = std::get<0>(fec); Vertex_handle u = std::get<0>(edge); Vertex_handle v = std::get<1>(edge); - Filtration_value filt = std::get<1>(fec); + Filtration filt = std::get<1>(fec); // Inserts the edge in the sparse matrix to update the graph (G_i) insert_new_edges(u, v, filt); diff --git a/src/Collapse/test/collapse_unit_test.cpp b/src/Collapse/test/collapse_unit_test.cpp index 38adfa8a..3a07e088 100644 --- a/src/Collapse/test/collapse_unit_test.cpp +++ b/src/Collapse/test/collapse_unit_test.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #define BOOST_TEST_DYN_LINK #define BOOST_TEST_MODULE "collapse" @@ -32,7 +33,7 @@ bool find_edge_in_list(const Filtered_edge& edge, const Filtered_sorted_edge_lis } return false; } - +/* void trace_and_check_collapse(const Filtered_sorted_edge_list& edges, const Filtered_sorted_edge_list& removed_edges) { std::cout << "BEFORE COLLAPSE - Total number of edges: " << edges.size() << std::endl; BOOST_CHECK(edges.size() > 0); @@ -68,70 +69,104 @@ void trace_and_check_collapse(const Filtered_sorted_edge_list& edges, const Filt } BOOST_AUTO_TEST_CASE(collapse) { - /* - 1 2 - o---o - | | - | | - | | - o---o - 0 3 - */ + // 1 2 + // o---o + // | | + // | | + // | | + // o---o + // 0 3 Filtered_sorted_edge_list edges {{1., 0, 1}, {1., 1, 2}, {1., 2, 3}, {1., 3, 0}}; trace_and_check_collapse(edges, {}); - /* - 1 2 - o---o - |\ /| - | x | - |/ \| - o---o - 0 3 - */ + // 1 2 + // o---o + // |\ /| + // | x | + // |/ \| + // o---o + // 0 3 edges.push_back({2., 0, 2}); edges.push_back({2., 1, 3}); trace_and_check_collapse(edges, {{2., 1, 3}}); - /* - 1 2 4 - o---o---o - |\ /| | - | x | | - |/ \| | - o---o---o - 0 3 5 - */ + // 1 2 4 + // o---o---o + // |\ /| | + // | x | | + // |/ \| | + // o---o---o + // 0 3 5 edges.push_back({3., 2, 4}); edges.push_back({3., 4, 5}); edges.push_back({3., 5, 3}); trace_and_check_collapse(edges, {{2., 1, 3}}); - /* - 1 2 4 - o---o---o - |\ /|\ /| - | x | x | - |/ \|/ \| - o---o---o - 0 3 5 - */ + // 1 2 4 + // o---o---o + // |\ /|\ /| + // | x | x | + // |/ \|/ \| + // o---o---o + // 0 3 5 edges.push_back({4., 2, 5}); edges.push_back({4., 4, 3}); trace_and_check_collapse(edges, {{2., 1, 3}, {4., 4, 3}}); - /* - 1 2 4 - o---o---o - |\ /|\ /| - | x | x | + [0,4] and [1,5] - |/ \|/ \| - o---o---o - 0 3 5 - */ + // 1 2 4 + // o---o---o + // |\ /|\ /| + // | x | x | + [0,4] and [1,5] + // |/ \|/ \| + // o---o---o + // 0 3 5 edges.push_back({5., 1, 5}); edges.push_back({5., 0, 4}); trace_and_check_collapse(edges, {{2., 1, 3}, {4., 4, 3}, {5., 0, 4}}); -} +}*/ + + +BOOST_AUTO_TEST_CASE(collapse_from_distance_matrix) { + // 1 2 + // o---o + // |\ /| + // | x | + // |/ \| + // o---o + // 0 3 + // Lower diagonal distance matrix + std::array, 4> distance_matrix = {{{0., 0., 0., 0.}, + {1., 0., 0., 0.}, + {2., 1., 0., 0.}, + {1., 2., 1., 0.} }}; + + std::cout << "COLLAPSE - keep edges: " << std::endl; + Flag_complex_sparse_matrix flag_complex_sparse_matrix(distance_matrix); + Filtered_sorted_edge_list collapse_edges; + flag_complex_sparse_matrix.filtered_edge_collapse( + [&collapse_edges](std::pair edge, Filtration_value filtration) { + std::cout << "f[" << std::get<0>(edge) << ", " << std::get<1>(edge) << "] = " << filtration << std::endl; + collapse_edges.push_back({filtration, std::get<0>(edge), std::get<1>(edge)}); + }); + std::cout << "AFTER COLLAPSE - Total number of edges: " << collapse_edges.size() << std::endl; + BOOST_CHECK(collapse_edges.size() == 5); + Filtered_sorted_edge_list edges {{1., 0, 1}, {1., 1, 2}, {1., 2, 3}, {1., 0, 3}, {2., 0, 2}, {2., 1, 3}}; + + for (auto edge_from_collapse : collapse_edges) { + std::cout << "f[" << std::get<1>(edge_from_collapse) << ", " << std::get<2>(edge_from_collapse) << "] = " + << std::get<0>(edge_from_collapse) << std::endl; + // Check each edge from collapse is in the input + BOOST_CHECK(find_edge_in_list(edge_from_collapse, edges)); + } + Filtered_sorted_edge_list removed_edges {{2., 1, 3}}; + + std::cout << "CHECK COLLAPSE - Total number of removed edges: " << removed_edges.size() << std::endl; + for (auto removed_edge : removed_edges) { + std::cout << "f[" << std::get<1>(removed_edge) << ", " << std::get<2>(removed_edge) << "] = " + << std::get<0>(removed_edge) << std::endl; + // Check each removed edge from collapse is in the input + BOOST_CHECK(!find_edge_in_list(removed_edge, collapse_edges)); + } +} \ No newline at end of file diff --git a/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp b/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp index f6926224..f4a460ab 100644 --- a/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp +++ b/src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp @@ -11,10 +11,8 @@ #include #include #include -#include -#include #include -#include +#include #include @@ -23,8 +21,8 @@ using Filtration_value = Simplex_tree::Filtration_value; using Vertex_handle = Simplex_tree::Vertex_handle; using Flag_complex_sparse_matrix = Gudhi::collapse::Flag_complex_sparse_matrix; +using Proximity_graph = Flag_complex_sparse_matrix::Proximity_graph; -using Rips_edge_list = Gudhi::rips_edge_list::Rips_edge_list; using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; using Distance_matrix = std::vector>; @@ -82,28 +80,22 @@ int main(int argc, char* argv[]) { program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence); Distance_matrix distances; - Distance_matrix sparse_distances; distances = Gudhi::read_lower_triangular_matrix_from_csv_file(csv_matrix_file); std::cout << "Read the distance matrix succesfully, of size: " << distances.size() << std::endl; - Flag_complex_sparse_matrix::Filtered_sorted_edge_list edge_t; - - Rips_edge_list Rips_edge_list_from_file(distances, threshold); - Rips_edge_list_from_file.create_edges(edge_t); - - if (edge_t.size() <= 0) { - std::cerr << "Total number of egdes are zero." << std::endl; - exit(-1); - } - - std::cout << "Total number of edges before collapse are: " << edge_t.size() << std::endl; + Proximity_graph proximity_graph = Gudhi::compute_proximity_graph(boost::irange((size_t)0, + distances.size()), + threshold, + [&distances](size_t i, size_t j) { + return distances[j][i]; + }); // Now we will perform filtered edge collapse to sparsify the edge list edge_t. - Flag_complex_sparse_matrix mat_filt_edge_coll(edge_t); + Flag_complex_sparse_matrix flag_complex(proximity_graph); Simplex_tree stree; - mat_filt_edge_coll.filtered_edge_collapse( + flag_complex.filtered_edge_collapse( [&stree](std::vector edge, Filtration_value filtration) { // insert the 2 vertices with a 0. filtration value just like a Rips stree.insert_simplex({edge[0]}, 0.); diff --git a/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp b/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp index e322d3cd..b9130d4c 100644 --- a/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp +++ b/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp @@ -12,11 +12,9 @@ #include #include #include -#include #include #include -#include #include #include // for std::pair @@ -31,15 +29,12 @@ using Point = std::vector; using Vector_of_points = std::vector; using Flag_complex_sparse_matrix = Gudhi::collapse::Flag_complex_sparse_matrix; +using Proximity_graph = Flag_complex_sparse_matrix::Proximity_graph; using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; using Distance_matrix = std::vector>; -using Adjacency_list = boost::adjacency_list, - boost::property>; - void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence); @@ -76,9 +71,9 @@ int main(int argc, char* argv[]) { std::cout << "Successfully read " << point_vector.size() << " point_vector.\n"; std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n"; - Adjacency_list proximity_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), - threshold, - Gudhi::Euclidean_distance()); + Proximity_graph proximity_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), + threshold, + Gudhi::Euclidean_distance()); if (num_edges(proximity_graph) <= 0) { std::cerr << "Total number of egdes are zero." << std::endl; @@ -89,7 +84,7 @@ int main(int argc, char* argv[]) { Simplex_tree stree; mat_filt_edge_coll.filtered_edge_collapse( - [&stree](std::vector edge, Filtration_value filtration) { + [&stree](const std::vector& edge, Filtration_value filtration) { // insert the 2 vertices with a 0. filtration value just like a Rips stree.insert_simplex({edge[0]}, 0.); stree.insert_simplex({edge[1]}, 0.); -- cgit v1.2.3