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