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>
23 #include <tbb/parallel_sort.h>
29 #include <unordered_map>
30 #include <unordered_set>
35 #include <type_traits>
52 template<
typename Vertex,
typename Filtration>
53 class Flag_complex_edge_collapser {
62 using IVertex = std::size_t;
63 using Edge_index = std::size_t;
64 using IEdge = std::pair<IVertex, IVertex>;
68 using Sparse_vector = Eigen::SparseVector<Edge_index>;
69 using Sparse_row_matrix = std::vector<Sparse_vector>;
74 class iterator :
public boost::iterator_facade<iterator,
76 std::input_iterator_tag,
80 iterator():ptr(nullptr){}
81 iterator(Neighbours
const*p):ptr(p){find_valid();}
83 friend class boost::iterator_core_access;
92 if(!it) { ptr=
nullptr;
break; }
93 if(IVertex(it.index()) == ptr->u) {
97 Edge_index e = it.value();
98 if(e <= ptr->ec->current_backward || ptr->ec->critical_edge_indicator_[e])
break;
101 bool equal(iterator
const& other)
const {
return ptr == other.ptr; }
102 IVertex dereference()
const {
return ptr->it.index(); }
104 typedef iterator const_iterator;
105 mutable typename Sparse_vector::InnerIterator it;
106 Flag_complex_edge_collapser
const*ec;
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){}
114 using IVertex_vector = std::vector<IVertex>;
118 using Filtered_edge = std::tuple<Vertex_handle, Vertex_handle, Filtration_value>;
122 std::vector<Vertex_handle> row_to_vertex_;
126 Edge_index current_backward = -1;
129 std::unordered_map<IEdge, Edge_index, boost::hash<IEdge>> iedge_to_index_map_;
132 std::vector<bool> critical_edge_indicator_;
135 std::unordered_map<Vertex_handle, IVertex> vertex_to_row_;
139 Sparse_row_matrix sparse_row_adjacency_matrix_;
142 std::vector<Filtered_edge> f_edge_vector_;
147 const IVertex rw_u = vertex_to_row_.at(u);
148 const IVertex rw_v = vertex_to_row_.at(v);
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);
154 std::cout <<
"And its common neighbours are." << std::endl;
155 for (
auto neighbour : common_neighbours) {
156 std::cout << row_to_vertex_[neighbour] <<
", " ;
158 std::cout<< std::endl;
159 #endif // DEBUG_TRACES
160 if (common_neighbours.size() == 1)
163 for (
auto rw_c : common_neighbours) {
164 auto neighbours_c = neighbours_row_index<true>(rw_c);
166 if (std::includes(neighbours_c.begin(), neighbours_c.end(),
167 common_neighbours.begin(), common_neighbours.end()))
174 std::set<Edge_index> three_clique_indices(Edge_index crit) {
175 std::set<Edge_index> edge_indices;
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];
187 IVertex_vector common_neighbours = open_common_neighbours_row_index(rw_u, rw_v);
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]);
199 template<
typename FilteredEdgeOutput>
200 void set_edge_critical(Edge_index indx,
Filtration_value filt, FilteredEdgeOutput filtered_edge_output) {
202 std::cout <<
"The curent index with filtration value " << indx <<
", " << filt <<
" is primary critical" <<
204 #endif // DEBUG_TRACES
205 std::set<Edge_index> effected_indices = three_clique_indices(indx);
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]);
212 if (!critical_edge_indicator_[current_backward]) {
213 if (!edge_is_dominated(u, v)) {
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)
222 effected_indices.emplace(inr_idx);
225 std::cout <<
"The following edge is critical with filt value: {" << u <<
"," << v <<
"}; "
226 << filt << std::endl;
227 #endif // DEBUG_TRACES
232 current_backward = -1;
236 template<
bool closed>
237 auto neighbours_row_index(IVertex rw_u)
const
239 return Neighbours<closed>(
this, rw_u);
243 IVertex_vector open_common_neighbours_row_index(IVertex rw_u, IVertex rw_v)
const
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));
256 auto n = row_to_vertex_.size();
257 auto result = vertex_to_row_.emplace(vertex, n);
261 sparse_row_adjacency_matrix_.emplace_back((std::numeric_limits<Eigen::Index>::max)());
263 sparse_row_adjacency_matrix_[n].insert(n) = -1;
265 row_to_vertex_.push_back(vertex);
267 return result.first->second;
274 GUDHI_CHECK((u != v), std::invalid_argument(
"Flag_complex_edge_collapser::insert_new_edge with u == v"));
276 IVertex rw_u = insert_vertex(u);
277 IVertex rw_v = insert_vertex(v);
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);
295 template<
typename FilteredEdgeRange>
296 Flag_complex_edge_collapser(
const FilteredEdgeRange& edges)
297 : f_edge_vector_(std::begin(edges), std::end(edges)) { }
305 template<
typename FilteredEdgeOutput>
306 void process_edges(FilteredEdgeOutput filtered_edge_output) {
308 auto sort_by_filtration = [](
const Filtered_edge& edge_a,
const Filtered_edge& edge_b) ->
bool
310 return (std::get<2>(edge_a) < std::get<2>(edge_b));
314 tbb::parallel_sort(f_edge_vector_.begin(), f_edge_vector_.end(), sort_by_filtration);
316 std::sort(f_edge_vector_.begin(), f_edge_vector_.end(), sort_by_filtration);
319 for (Edge_index endIdx = 0; endIdx < f_edge_vector_.size(); endIdx++) {
320 Filtered_edge fec = f_edge_vector_[endIdx];
326 IEdge ie = insert_new_edge(u, v, endIdx);
328 iedge_to_index_map_.emplace(ie, endIdx);
329 critical_edge_indicator_.push_back(
false);
331 if (!edge_is_dominated(u, v)) {
332 critical_edge_indicator_[endIdx] =
true;
333 filtered_edge_output(u, v, filt);
335 set_edge_critical(endIdx, filt, filtered_edge_output);
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(
368 remaining_edges.emplace_back(u, v, filtration);
371 return remaining_edges;
378 #endif // FLAG_COMPLEX_EDGE_COLLAPSER_H_