Flag_complex_edge_collapser.h
1/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2 * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3 * Author(s): Siddharth Pritam
4 *
5 * Copyright (C) 2020 Inria
6 *
7 * Modification(s):
8 * - 2020/03 Vincent Rouvreau: integration to the gudhi library
9 * - YYYY/MM Author: Description of the modification
10 */
11
12#ifndef FLAG_COMPLEX_EDGE_COLLAPSER_H_
13#define FLAG_COMPLEX_EDGE_COLLAPSER_H_
14
15#include <gudhi/Debug_utils.h>
16
17#include <boost/functional/hash.hpp>
18#include <boost/iterator/iterator_facade.hpp>
19
20#include <Eigen/Sparse>
21#include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST
22
23#ifdef GUDHI_USE_TBB
24#include <tbb/parallel_sort.h>
25#endif
26
27#include <iostream>
28#include <utility> // for std::pair
29#include <vector>
30#include <unordered_map>
31#include <unordered_set>
32#include <set>
33#include <tuple> // for std::tie
34#include <algorithm> // for std::includes
35#include <iterator> // for std::inserter
36#include <type_traits> // for std::decay
37
38// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
39#if !EIGEN_VERSION_AT_LEAST(3,1,0)
40# error Edge Collapse is only available for Eigen3 >= 3.1.0
41#endif
42
43namespace Gudhi {
44
45namespace collapse {
46
58template<typename Vertex, typename Filtration>
59class Flag_complex_edge_collapser {
60 public:
62 using Vertex_handle = Vertex;
64 using Filtration_value = Filtration;
65
66 private:
67 // internal numbering of vertices and edges
68 using IVertex = std::size_t;
69 using Edge_index = std::size_t;
70 using IEdge = std::pair<IVertex, IVertex>;
71
72 // The sparse matrix data type
73 // (Eigen::SparseMatrix<Edge_index, Eigen::RowMajor> has slow insertions)
74 using Sparse_vector = Eigen::SparseVector<Edge_index>;
75 using Sparse_row_matrix = std::vector<Sparse_vector>;
76
77 // Range of neighbors of a vertex
78 template<bool closed>
79 struct Neighbours {
80 class iterator : public boost::iterator_facade<iterator,
81 IVertex, /* value_type */
82 std::input_iterator_tag, // or boost::single_pass_traversal_tag
83 IVertex /* reference */ >
84 {
85 public:
86 iterator():ptr(nullptr){}
87 iterator(Neighbours const*p):ptr(p){find_valid();}
88 private:
89 friend class boost::iterator_core_access;
90 Neighbours const*ptr;
91 void increment(){
92 ++ptr->it;
93 find_valid();
94 }
95 void find_valid(){
96 auto& it = ptr->it;
97 do {
98 if(!it) { ptr=nullptr; break; }
99 if(IVertex(it.index()) == ptr->u) {
100 if(closed) break;
101 else continue;
102 }
103 Edge_index e = it.value();
104 if(e <= ptr->ec->current_backward || ptr->ec->critical_edge_indicator_[e]) break;
105 } while(++it, true);
106 }
107 bool equal(iterator const& other) const { return ptr == other.ptr; }
108 IVertex dereference() const { return ptr->it.index(); }
109 };
110 typedef iterator const_iterator;
111 mutable typename Sparse_vector::InnerIterator it;
112 Flag_complex_edge_collapser const*ec;
113 IVertex u;
114 iterator begin() const { return this; }
115 iterator end() const { return {}; }
116 explicit Neighbours(Flag_complex_edge_collapser const*p,IVertex u):it(p->sparse_row_adjacency_matrix_[u]),ec(p),u(u){}
117 };
118
119 // A range of row indices
120 using IVertex_vector = std::vector<IVertex>;
121
122 public:
124 using Filtered_edge = std::tuple<Vertex_handle, Vertex_handle, Filtration_value>;
125
126 private:
127 // Map from row index to its vertex handle
128 std::vector<Vertex_handle> row_to_vertex_;
129
130 // Index of the current edge in the backwards walk. Edges <= current_backward are part of the temporary graph,
131 // while edges > current_backward are removed unless critical_edge_indicator_.
132 Edge_index current_backward = -1;
133
134 // Map from IEdge to its index
135 std::unordered_map<IEdge, Edge_index, boost::hash<IEdge>> iedge_to_index_map_;
136
137 // Boolean vector to indicate if the edge is critical.
138 std::vector<bool> critical_edge_indicator_;
139
140 // Map from vertex handle to its row index
141 std::unordered_map<Vertex_handle, IVertex> vertex_to_row_;
142
143 // Stores the Sparse matrix of Filtration values representing the original graph.
144 // The matrix rows and columns are indexed by IVertex.
145 Sparse_row_matrix sparse_row_adjacency_matrix_;
146
147 // The input, a vector of filtered edges.
148 std::vector<Filtered_edge> f_edge_vector_;
149
150 // Edge is the actual edge (u,v), with Vertex_handle u and v, not IVertex.
151 bool edge_is_dominated(Vertex_handle u, Vertex_handle v) const
152 {
153 const IVertex rw_u = vertex_to_row_.at(u);
154 const IVertex rw_v = vertex_to_row_.at(v);
155#ifdef DEBUG_TRACES
156 std::cout << "The edge {" << u << ", " << v << "} is going for domination check." << std::endl;
157#endif // DEBUG_TRACES
158 auto common_neighbours = open_common_neighbours_row_index(rw_u, rw_v);
159#ifdef DEBUG_TRACES
160 std::cout << "And its common neighbours are." << std::endl;
161 for (auto neighbour : common_neighbours) {
162 std::cout << row_to_vertex_[neighbour] << ", " ;
163 }
164 std::cout<< std::endl;
165#endif // DEBUG_TRACES
166 if (common_neighbours.size() == 1)
167 return true;
168 else
169 for (auto rw_c : common_neighbours) {
170 auto neighbours_c = neighbours_row_index<true>(rw_c);
171 // If neighbours_c contains the common neighbours.
172 if (std::includes(neighbours_c.begin(), neighbours_c.end(),
173 common_neighbours.begin(), common_neighbours.end()))
174 return true;
175 }
176 return false;
177 }
178
179 // Returns the edges connecting u and v (extremities of crit) to their common neighbors (not themselves)
180 std::set<Edge_index> three_clique_indices(Edge_index crit) {
181 std::set<Edge_index> edge_indices;
182
183 Vertex_handle u = std::get<0>(f_edge_vector_[crit]);
184 Vertex_handle v = std::get<1>(f_edge_vector_[crit]);
185
186#ifdef DEBUG_TRACES
187 std::cout << "The current critical edge to re-check criticality with filt value is : f {" << u << "," << v
188 << "} = " << std::get<2>(f_edge_vector_[crit]) << std::endl;
189#endif // DEBUG_TRACES
190 auto rw_u = vertex_to_row_[u];
191 auto rw_v = vertex_to_row_[v];
192
193 IVertex_vector common_neighbours = open_common_neighbours_row_index(rw_u, rw_v);
194
195 for (auto rw_c : common_neighbours) {
196 IEdge e_with_new_nbhr_v = std::minmax(rw_u, rw_c);
197 IEdge e_with_new_nbhr_u = std::minmax(rw_v, rw_c);
198 edge_indices.emplace(iedge_to_index_map_[e_with_new_nbhr_v]);
199 edge_indices.emplace(iedge_to_index_map_[e_with_new_nbhr_u]);
200 }
201 return edge_indices;
202 }
203
204 // Detect and set all edges that are becoming critical
205 template<typename FilteredEdgeOutput>
206 void set_edge_critical(Edge_index indx, Filtration_value filt, FilteredEdgeOutput filtered_edge_output) {
207#ifdef DEBUG_TRACES
208 std::cout << "The curent index with filtration value " << indx << ", " << filt << " is primary critical" <<
209 std::endl;
210#endif // DEBUG_TRACES
211 std::set<Edge_index> effected_indices = three_clique_indices(indx);
212 // Cannot use boost::adaptors::reverse in such dynamic cases apparently
213 for (auto it = effected_indices.rbegin(); it != effected_indices.rend(); ++it) {
214 current_backward = *it;
215 Vertex_handle u = std::get<0>(f_edge_vector_[current_backward]);
216 Vertex_handle v = std::get<1>(f_edge_vector_[current_backward]);
217 // If current_backward is not critical so it should be processed, otherwise it stays in the graph
218 if (!critical_edge_indicator_[current_backward]) {
219 if (!edge_is_dominated(u, v)) {
220#ifdef DEBUG_TRACES
221 std::cout << "The curent index became critical " << current_backward << std::endl;
222#endif // DEBUG_TRACES
223 critical_edge_indicator_[current_backward] = true;
224 filtered_edge_output(u, v, filt);
225 std::set<Edge_index> inner_effected_indcs = three_clique_indices(current_backward);
226 for (auto inr_idx : inner_effected_indcs) {
227 if(inr_idx < current_backward) // && !critical_edge_indicator_[inr_idx]
228 effected_indices.emplace(inr_idx);
229 }
230#ifdef DEBUG_TRACES
231 std::cout << "The following edge is critical with filt value: {" << u << "," << v << "}; "
232 << filt << std::endl;
233#endif // DEBUG_TRACES
234 }
235 }
236 }
237 // Clear the implicit "removed from graph" data structure
238 current_backward = -1;
239 }
240
241 // Returns list of neighbors of a particular vertex.
242 template<bool closed>
243 auto neighbours_row_index(IVertex rw_u) const
244 {
245 return Neighbours<closed>(this, rw_u);
246 }
247
248 // Returns the list of open neighbours of the edge :{u,v}.
249 IVertex_vector open_common_neighbours_row_index(IVertex rw_u, IVertex rw_v) const
250 {
251 auto non_zero_indices_u = neighbours_row_index<false>(rw_u);
252 auto non_zero_indices_v = neighbours_row_index<false>(rw_v);
253 IVertex_vector common;
254 std::set_intersection(non_zero_indices_u.begin(), non_zero_indices_u.end(), non_zero_indices_v.begin(),
255 non_zero_indices_v.end(), std::back_inserter(common));
256
257 return common;
258 }
259
260 // Insert a vertex in the data structure
261 IVertex insert_vertex(Vertex_handle vertex) {
262 auto n = row_to_vertex_.size();
263 auto result = vertex_to_row_.emplace(vertex, n);
264 // If it was not already inserted - Value won't be updated by emplace if it is already present
265 if (result.second) {
266 // Expand the matrix. The size of rows is irrelevant.
267 sparse_row_adjacency_matrix_.emplace_back((std::numeric_limits<Eigen::Index>::max)());
268 // Initializing the diagonal element of the adjency matrix corresponding to rw_b.
269 sparse_row_adjacency_matrix_[n].insert(n) = -1; // not an edge
270 // Must be done after reading its size()
271 row_to_vertex_.push_back(vertex);
272 }
273 return result.first->second;
274 }
275
276 // Insert an edge in the data structure
277 // @exception std::invalid_argument In debug mode, if u == v
278 IEdge insert_new_edge(Vertex_handle u, Vertex_handle v, Edge_index idx)
279 {
280 GUDHI_CHECK((u != v), std::invalid_argument("Flag_complex_edge_collapser::insert_new_edge with u == v"));
281 // The edge must not be added before, it should be a new edge.
282 IVertex rw_u = insert_vertex(u);
283 IVertex rw_v = insert_vertex(v);
284#ifdef DEBUG_TRACES
285 std::cout << "Inserting the edge " << u <<", " << v << std::endl;
286#endif // DEBUG_TRACES
287 sparse_row_adjacency_matrix_[rw_u].insert(rw_v) = idx;
288 sparse_row_adjacency_matrix_[rw_v].insert(rw_u) = idx;
289 return std::minmax(rw_u, rw_v);
290 }
291
292 public:
301 template<typename FilteredEdgeRange>
302 Flag_complex_edge_collapser(const FilteredEdgeRange& edges)
303 : f_edge_vector_(std::begin(edges), std::end(edges)) { }
304
311 template<typename FilteredEdgeOutput>
312 void process_edges(FilteredEdgeOutput filtered_edge_output) {
313 // Sort edges
314 auto sort_by_filtration = [](const Filtered_edge& edge_a, const Filtered_edge& edge_b) -> bool
315 {
316 return (std::get<2>(edge_a) < std::get<2>(edge_b));
317 };
318
319#ifdef GUDHI_USE_TBB
320 tbb::parallel_sort(f_edge_vector_.begin(), f_edge_vector_.end(), sort_by_filtration);
321#else
322 std::sort(f_edge_vector_.begin(), f_edge_vector_.end(), sort_by_filtration);
323#endif
324
325 for (Edge_index endIdx = 0; endIdx < f_edge_vector_.size(); endIdx++) {
326 Filtered_edge fec = f_edge_vector_[endIdx];
327 Vertex_handle u = std::get<0>(fec);
328 Vertex_handle v = std::get<1>(fec);
329 Filtration_value filt = std::get<2>(fec);
330
331 // Inserts the edge in the sparse matrix to update the graph (G_i)
332 IEdge ie = insert_new_edge(u, v, endIdx);
333
334 iedge_to_index_map_.emplace(ie, endIdx);
335 critical_edge_indicator_.push_back(false);
336
337 if (!edge_is_dominated(u, v)) {
338 critical_edge_indicator_[endIdx] = true;
339 filtered_edge_output(u, v, filt);
340 if (endIdx > 1)
341 set_edge_critical(endIdx, filt, filtered_edge_output);
342 }
343 }
344 }
345
346};
347
363template<class FilteredEdgeRange> auto flag_complex_collapse_edges(const FilteredEdgeRange& edges) {
364 auto first_edge_itr = std::begin(edges);
365 using Vertex_handle = std::decay_t<decltype(std::get<0>(*first_edge_itr))>;
366 using Filtration_value = std::decay_t<decltype(std::get<2>(*first_edge_itr))>;
367 using Edge_collapser = Flag_complex_edge_collapser<Vertex_handle, Filtration_value>;
368 std::vector<typename Edge_collapser::Filtered_edge> remaining_edges;
369 if (first_edge_itr != std::end(edges)) {
370 Edge_collapser edge_collapser(edges);
371 edge_collapser.process_edges(
372 [&remaining_edges](Vertex_handle u, Vertex_handle v, Filtration_value filtration) {
373 // insert the edge
374 remaining_edges.emplace_back(u, v, filtration);
375 });
376 }
377 return remaining_edges;
378}
379
380} // namespace collapse
381
382} // namespace Gudhi
383
384#endif // FLAG_COMPLEX_EDGE_COLLAPSER_H_
auto flag_complex_collapse_edges(const FilteredEdgeRange &edges)
Implicitly constructs a flag complex from edges as an input, collapses edges while preserving the per...
Definition: Flag_complex_edge_collapser.h:363
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20
Handle type for the vertices of a cell complex.
Definition: VertexHandle.h:15
GUDHIdev  Version 3.5.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Tue Aug 16 2022 14:01:50 for GUDHIdev by Doxygen 1.9.4