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 
43 namespace Gudhi {
44 
45 namespace collapse {
46 
58 template<typename Vertex, typename Filtration>
59 class 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 
363 template<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_
Gudhi::collapse::flag_complex_collapse_edges
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
VertexHandle
Handle type for the vertices of a cell complex.
Definition: VertexHandle.h:15
FiltrationValue
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20
GUDHI  Version 3.4.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Tue Dec 29 2020 10:00:12 for GUDHI by Doxygen 1.8.20