12 #ifndef FLAG_COMPLEX_EDGE_COLLAPSER_H_
13 #define FLAG_COMPLEX_EDGE_COLLAPSER_H_
15 #include <gudhi/Debug_utils.h>
17 #include <boost/functional/hash.hpp>
18 #include <boost/iterator/iterator_facade.hpp>
20 #include <Eigen/Sparse>
21 #include <Eigen/src/Core/util/Macros.h>
24 #include <tbb/parallel_sort.h>
30 #include <unordered_map>
31 #include <unordered_set>
36 #include <type_traits>
39 #if !EIGEN_VERSION_AT_LEAST(3,1,0)
40 # error Edge Collapse is only available for Eigen3 >= 3.1.0
58 template<
typename Vertex,
typename Filtration>
59 class Flag_complex_edge_collapser {
68 using IVertex = std::size_t;
69 using Edge_index = std::size_t;
70 using IEdge = std::pair<IVertex, IVertex>;
74 using Sparse_vector = Eigen::SparseVector<Edge_index>;
75 using Sparse_row_matrix = std::vector<Sparse_vector>;
80 class iterator :
public boost::iterator_facade<iterator,
82 std::input_iterator_tag,
86 iterator():ptr(nullptr){}
87 iterator(Neighbours
const*p):ptr(p){find_valid();}
89 friend class boost::iterator_core_access;
98 if(!it) { ptr=
nullptr;
break; }
99 if(IVertex(it.index()) == ptr->u) {
103 Edge_index e = it.value();
104 if(e <= ptr->ec->current_backward || ptr->ec->critical_edge_indicator_[e])
break;
107 bool equal(iterator
const& other)
const {
return ptr == other.ptr; }
108 IVertex dereference()
const {
return ptr->it.index(); }
110 typedef iterator const_iterator;
111 mutable typename Sparse_vector::InnerIterator it;
112 Flag_complex_edge_collapser
const*ec;
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){}
120 using IVertex_vector = std::vector<IVertex>;
124 using Filtered_edge = std::tuple<Vertex_handle, Vertex_handle, Filtration_value>;
128 std::vector<Vertex_handle> row_to_vertex_;
132 Edge_index current_backward = -1;
135 std::unordered_map<IEdge, Edge_index, boost::hash<IEdge>> iedge_to_index_map_;
138 std::vector<bool> critical_edge_indicator_;
141 std::unordered_map<Vertex_handle, IVertex> vertex_to_row_;
145 Sparse_row_matrix sparse_row_adjacency_matrix_;
148 std::vector<Filtered_edge> f_edge_vector_;
153 const IVertex rw_u = vertex_to_row_.at(u);
154 const IVertex rw_v = vertex_to_row_.at(v);
156 std::cout <<
"The edge {" << u <<
", " << v <<
"} is going for domination check." << std::endl;
158 auto common_neighbours = open_common_neighbours_row_index(rw_u, rw_v);
160 std::cout <<
"And its common neighbours are." << std::endl;
161 for (
auto neighbour : common_neighbours) {
162 std::cout << row_to_vertex_[neighbour] <<
", " ;
164 std::cout<< std::endl;
166 if (common_neighbours.size() == 1)
169 for (
auto rw_c : common_neighbours) {
170 auto neighbours_c = neighbours_row_index<true>(rw_c);
172 if (std::includes(neighbours_c.begin(), neighbours_c.end(),
173 common_neighbours.begin(), common_neighbours.end()))
180 std::set<Edge_index> three_clique_indices(Edge_index crit) {
181 std::set<Edge_index> edge_indices;
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;
190 auto rw_u = vertex_to_row_[u];
191 auto rw_v = vertex_to_row_[v];
193 IVertex_vector common_neighbours = open_common_neighbours_row_index(rw_u, rw_v);
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]);
205 template<
typename FilteredEdgeOutput>
206 void set_edge_critical(Edge_index indx,
Filtration_value filt, FilteredEdgeOutput filtered_edge_output) {
208 std::cout <<
"The curent index with filtration value " << indx <<
", " << filt <<
" is primary critical" <<
211 std::set<Edge_index> effected_indices = three_clique_indices(indx);
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]);
218 if (!critical_edge_indicator_[current_backward]) {
219 if (!edge_is_dominated(u, v)) {
221 std::cout <<
"The curent index became critical " << current_backward << std::endl;
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)
228 effected_indices.emplace(inr_idx);
231 std::cout <<
"The following edge is critical with filt value: {" << u <<
"," << v <<
"}; "
232 << filt << std::endl;
238 current_backward = -1;
242 template<
bool closed>
243 auto neighbours_row_index(IVertex rw_u)
const
245 return Neighbours<closed>(
this, rw_u);
249 IVertex_vector open_common_neighbours_row_index(IVertex rw_u, IVertex rw_v)
const
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));
262 auto n = row_to_vertex_.size();
263 auto result = vertex_to_row_.emplace(vertex, n);
267 sparse_row_adjacency_matrix_.emplace_back((std::numeric_limits<Eigen::Index>::max)());
269 sparse_row_adjacency_matrix_[n].insert(n) = -1;
271 row_to_vertex_.push_back(vertex);
273 return result.first->second;
280 GUDHI_CHECK((u != v), std::invalid_argument(
"Flag_complex_edge_collapser::insert_new_edge with u == v"));
282 IVertex rw_u = insert_vertex(u);
283 IVertex rw_v = insert_vertex(v);
285 std::cout <<
"Inserting the edge " << u <<
", " << v << std::endl;
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);
301 template<
typename FilteredEdgeRange>
302 Flag_complex_edge_collapser(
const FilteredEdgeRange& edges)
303 : f_edge_vector_(std::begin(edges), std::end(edges)) { }
311 template<
typename FilteredEdgeOutput>
312 void process_edges(FilteredEdgeOutput filtered_edge_output) {
314 auto sort_by_filtration = [](
const Filtered_edge& edge_a,
const Filtered_edge& edge_b) ->
bool
316 return (std::get<2>(edge_a) < std::get<2>(edge_b));
320 tbb::parallel_sort(f_edge_vector_.begin(), f_edge_vector_.end(), sort_by_filtration);
322 std::sort(f_edge_vector_.begin(), f_edge_vector_.end(), sort_by_filtration);
325 for (Edge_index endIdx = 0; endIdx < f_edge_vector_.size(); endIdx++) {
326 Filtered_edge fec = f_edge_vector_[endIdx];
332 IEdge ie = insert_new_edge(u, v, endIdx);
334 iedge_to_index_map_.emplace(ie, endIdx);
335 critical_edge_indicator_.push_back(
false);
337 if (!edge_is_dominated(u, v)) {
338 critical_edge_indicator_[endIdx] =
true;
339 filtered_edge_output(u, v, filt);
341 set_edge_critical(endIdx, filt, filtered_edge_output);
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(
374 remaining_edges.emplace_back(u, v, filtration);
377 return remaining_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
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