summaryrefslogtreecommitdiff
path: root/src/Collapse
diff options
context:
space:
mode:
authorROUVREAU Vincent <vincent.rouvreau@inria.fr>2020-04-12 10:33:38 +0200
committerROUVREAU Vincent <vincent.rouvreau@inria.fr>2020-04-12 10:33:38 +0200
commit1ce5d0d19e13a14e8a67442aec7bc40eae68dc8e (patch)
tree4e4c148cfd7935e8a686a7cafaaaacba2db432b9 /src/Collapse
parentf000725296c1962155e2ec331a3db6244d7c9f9e (diff)
Modify interface and utils
Diffstat (limited to 'src/Collapse')
-rw-r--r--src/Collapse/include/gudhi/Flag_complex_sparse_matrix.h86
-rw-r--r--src/Collapse/test/collapse_unit_test.cpp129
-rw-r--r--src/Collapse/utilities/distance_matrix_edge_collapse_rips_persistence.cpp28
-rw-r--r--src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp15
4 files changed, 148 insertions, 110 deletions
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 <gudhi/graph_simplicial_complex.h>
#include <boost/functional/hash.hpp>
+#include <boost/graph/adjacency_list.hpp>
#include <Eigen/Sparse>
@@ -32,6 +33,7 @@
#include <list>
#include <algorithm>
+
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<typename Vertex_handle, typename Filtration_value>
+template<typename Vertex_handle, typename Filtration>
class Flag_complex_sparse_matrix {
private:
using Edge = std::pair<Vertex_handle, Vertex_handle>; // 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<Edge, Filtration_value>;
+ using Filtered_edge = std::pair<Edge, Filtration>;
using Row_index = std::size_t;
using Map_vertex_to_index = std::unordered_map<Vertex_handle, Row_index>;
- using Sparse_row_matrix = Eigen::SparseMatrix<Filtration_value, Eigen::RowMajor>;
+ using Sparse_row_matrix = Eigen::SparseMatrix<Filtration, Eigen::RowMajor>;
using Row_indices_vector = std::vector<Row_index>;
public:
- using Filtered_sorted_edge_list = std::vector<std::tuple<Filtration_value, Vertex_handle, Vertex_handle>>;
+ using Filtered_sorted_edge_list = std::vector<std::tuple<Filtration, Vertex_handle, Vertex_handle>>;
+ using Filtration_value = Filtration;
+ using Proximity_graph = Gudhi::Proximity_graph<Flag_complex_sparse_matrix>;
private:
std::unordered_map<Row_index, Vertex_handle> 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<bool> dominated_edge_indicator_; // domination indicator
- //! Stores the Map between vertices_<B>row_to_vertex_ and row indices <B>row_to_vertex_ -> row-index</B>.
+ //! Stores the Map between vertices<B>row_to_vertex_ and row indices <B>row_to_vertex_ -> row-index</B>.
/*!
- So, if the original simplex tree had vertices_ 0,1,4,5 <br>
+ So, if the original simplex tree had vertices 0,1,4,5 <br>
<B>row_to_vertex_</B> would store : <br>
\verbatim
Values = | 0 | 1 | 4 | 5 |
@@ -106,10 +110,10 @@ class Flag_complex_sparse_matrix {
*/
std::unordered_map<Vertex_handle, Row_index> 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<Filtration_value, Eigen::RowMajor> ;
+ Sparse_row_matrix = Eigen::SparseMatrix<Filtration, Eigen::RowMajor> ;
\endcode
;
*/
@@ -120,7 +124,7 @@ class Flag_complex_sparse_matrix {
//! Stores <I>true</I> for dominated rows and <I>false</I> for undominated rows.
/*!
Initialised to a vector of length equal to the value of the variable <B>rows</B> with all <I>false</I> values.
- Subsequent removal of dominated vertices_ is reflected by concerned entries changing to <I>true</I> in this vector.
+ Subsequent removal of dominated vertices is reflected by concerned entries changing to <I>true</I> in this vector.
*/
std::vector<bool> domination_indicator_; //(domination indicator)
@@ -128,9 +132,9 @@ class Flag_complex_sparse_matrix {
std::vector<Filtered_edge> 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<typename FilteredEdgeInsertion>
- 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 {
<B>domination_indicator_</B> are initialised by init() function which is
called at the begining of this. <br>
*/
- Flag_complex_sparse_matrix(const Filtered_sorted_edge_list& edge_t)
+ template<typename DistanceMatrix>
+ 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<typename OneSkeletonGraph>
- 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<typename FilteredEdgeInsertion>
+ 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<typename FilteredEdgeInsertion>
- 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 <iostream>
#include <tuple>
#include <vector>
+#include <array>
#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<std::array<double, 4>, 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<Vertex_handle, Vertex_handle> 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 <gudhi/Flag_complex_sparse_matrix.h>
#include <gudhi/Simplex_tree.h>
#include <gudhi/Persistent_cohomology.h>
-#include <gudhi/Rips_edge_list.h>
-#include <gudhi/distance_functions.h>
#include <gudhi/reader_utils.h>
-#include <gudhi/Points_off_io.h>
+#include <gudhi/graph_simplicial_complex.h>
#include <boost/program_options.hpp>
@@ -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<Vertex_handle, Filtration_value>;
+using Proximity_graph = Flag_complex_sparse_matrix::Proximity_graph;
-using Rips_edge_list = Gudhi::rips_edge_list::Rips_edge_list<Filtration_value>;
using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>;
using Distance_matrix = std::vector<std::vector<Filtration_value>>;
@@ -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<Filtration_value>(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<Simplex_tree>(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<Vertex_handle> 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 <gudhi/Simplex_tree.h>
#include <gudhi/Persistent_cohomology.h>
#include <gudhi/distance_functions.h>
-#include <gudhi/reader_utils.h>
#include <gudhi/Points_off_io.h>
#include <gudhi/graph_simplicial_complex.h>
-#include <boost/graph/adjacency_list.hpp>
#include <boost/program_options.hpp>
#include<utility> // for std::pair
@@ -31,15 +29,12 @@ using Point = std::vector<Filtration_value>;
using Vector_of_points = std::vector<Point>;
using Flag_complex_sparse_matrix = Gudhi::collapse::Flag_complex_sparse_matrix<Vertex_handle, Filtration_value>;
+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<Simplex_tree, Field_Zp>;
using Distance_matrix = std::vector<std::vector<Filtration_value>>;
-using Adjacency_list = boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS,
- boost::property<Gudhi::vertex_filtration_t, double>,
- boost::property<Gudhi::edge_filtration_t, double>>;
-
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<Simplex_tree>(off_reader.get_point_cloud(),
- threshold,
- Gudhi::Euclidean_distance());
+ Proximity_graph proximity_graph = Gudhi::compute_proximity_graph<Simplex_tree>(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<Vertex_handle> edge, Filtration_value filtration) {
+ [&stree](const std::vector<Vertex_handle>& 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.);