From 2f7f48d2ca18356ba3430f6a30d2c8944c632be9 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Fri, 23 Jul 2021 00:02:05 +0200 Subject: Implement the backward algorithm --- .../include/gudhi/Flag_complex_edge_collapser.h | 247 ++++++++++++++++++++- 1 file changed, 246 insertions(+), 1 deletion(-) (limited to 'src/Collapse/include/gudhi/Flag_complex_edge_collapser.h') diff --git a/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h b/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h index 713c6608..3e0f2948 100644 --- a/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h +++ b/src/Collapse/include/gudhi/Flag_complex_edge_collapser.h @@ -1,6 +1,6 @@ /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. - * Author(s): Siddharth Pritam + * Author(s): Siddharth Pritam, Marc Glisse * * Copyright (C) 2020 Inria * @@ -16,6 +16,8 @@ #include #include +#include +#include #include #include // for EIGEN_VERSION_AT_LEAST @@ -345,6 +347,225 @@ class Flag_complex_edge_collapser { }; +template +struct Flag_complex_edge_collapser2 { + using Filtered_edge = std::tuple; + typedef std::pair Edge; + struct Cmpi { template bool operator()(T const&a, U const&b)const{return b Edge_group; // vector? on erase elements... + typedef std::map Events; // decreasing order + Events events; + typedef boost::container::flat_map Ngb_list; + //typedef std::unordered_map Neighbors; + typedef std::vector Neighbors; // we still need to find the edge in the group then... + Neighbors neighbors; + Vertex num_vertices; + + template + void read_edges(FilteredEdgeRange const&r){ + neighbors.resize(num_vertices); + std::vector neighbors2(num_vertices); + for(auto&&e : r){ + using std::get; + Vertex u = get<0>(e); + Vertex v = get<1>(e); + Filtration_value f = get<2>(e); + auto i = events.emplace_hint(events.end(), f, 0); + neighbors2[u].emplace_back(v, f); + neighbors2[v].emplace_back(u, f); + i->second.push_back(std::minmax({u, v})); + } + for(int i=0;i& e_ngb, std::vector>& e_ngb_later, Ngb_list const&nu, Ngb_list const&nv, Filtration_value f_event){ + auto ui = nu.begin(); + auto ue = nu.end(); + auto vi = nv.begin(); + auto ve = nv.end(); + assert(ui != ue && vi != ve); + while(ui != ue && vi != ve){ + Vertex w = ui->first; + if(w < vi->first) { ++ui; continue; } + if(w > vi->first) { ++vi; continue; } + Filtration_value f = std::max(ui->second, vi->second); + if(f > f_event) + e_ngb_later.emplace_back(f, w); + else + e_ngb.insert(e_ngb.end(), w); + ++ui; ++vi; + } + } + + // Test if the neighborhood of e is included in the closed neighborhood of c + template + bool is_dominated_by(Ngb const& e_ngb, Vertex c, Filtration_value f){ + Ngb_list const&nc = neighbors[c]; + // if few neighbors, use dichotomy? + // try a gallop strategy? an unordered_map? a bitset? +#if 0 + auto ci = nc.begin(); + auto ce = nc.end(); + for(auto v : e_ngb) { + if(v==c)continue; + ci = std::lower_bound(ci, ce, v, [](auto a, auto b){return a.first < b;}); + if(ci == nc.end() || ci->first != v || ci->second > f) return false; + } + return true; +#elif 0 + for(auto v : e_ngb) { + if(v==c)continue; + auto it = nc.find(v); + if(it == nc.end() || it->second > f) return false; + } + return true; +#else + auto ci = nc.begin(); + auto ce = nc.end(); + auto eni = e_ngb.begin(); + auto ene = e_ngb.end(); + assert(eni != ene); + assert(ci != ce); + if(*eni == c && ++eni == ene) return true; + for(;;){ + Vertex ve = *eni; + Vertex vc = ci->first; + while(ve > vc) { + if(++ci == ce) return false; + vc = ci->first; + } + if(ve < vc) return false; + // ve == vc + if(ci->second > f) return false; + if(++eni == ene)return true; + // Should we store a closed neighborhood of c (including c) so we can avoid testing for c at each iteration? + if(*eni == c && ++eni == ene)return true; + if(++ci == ce) return false; + } +#endif + } + + // delay = [](double d){return d*1.05;} + template + void process_edges(FilteredEdgeRange const& edges, Delay&& delay) { + { + Vertex maxi = 0, maxj = 0; + for(auto& fe : edges) { + Vertex i = std::get<0>(fe); + Vertex j = std::get<1>(fe); + if (i > maxi) maxi = i; + if (j > maxj) maxj = j; + } + num_vertices = std::max(maxi, maxj) + 1; + } + + read_edges(edges); + + typename Events::iterator evi; + boost::container::flat_set e_ngb; + e_ngb.reserve(num_vertices); + //std::multimap e_ngb_later; + std::vector> e_ngb_later; + // do not use reverse_iterator, it behaves badly with erase. + for(auto next_evi = events.begin(); (evi = next_evi) != events.end(); ){ + ++next_evi; // do it now, in case we remove evi + Edge_group& eg = evi->second; + //Filtration_value f_event = evi->first; + auto next_ei = eg.begin(); + auto ei = next_ei; + for(; (ei = next_ei) != eg.end();){ + next_ei = std::next(ei); // do it now, in case we remove ei + Vertex u = ei->first; + Vertex v = ei->second; + auto time = evi; + Filtration_value a_bit_later = delay(evi->first); + if(a_bit_later != evi->first) { +#if 0 + auto next_time = time; + while (next_time != events.begin() && (--next_time)->first <= a_bit_later) + time = next_time; +#else + // too bad there isn't a map::lower_bound_after + time = events.lower_bound(a_bit_later); +#endif + } + auto start_time = time; + e_ngb.clear(); + e_ngb_later.clear(); + common_neighbors(e_ngb, e_ngb_later, neighbors[u], neighbors[v], time->first); + auto cmp1=[](auto const&a, auto const&b){return a.first > b.first;}; + auto e_ngb_later_begin=e_ngb_later.begin(); + auto e_ngb_later_end=e_ngb_later.end(); + bool heapified = false; + + bool dead = false; + while(true) { + Vertex dominator = -1; + // special case for size 1 + // if(e_ngb.size()==1){dominator=*e_ngb.begin();}else + // TODO: try testing dominators in order of filtration value + for(auto c : e_ngb){ + if(is_dominated_by(e_ngb, c, time->first)){ + dominator = c; + break; + } + } + if(dominator==-1) break; + bool still_dominated; + do { + if(e_ngb_later_begin == e_ngb_later_end) { + dead = true; goto end_move; + } + if(!heapified) { + // Eagerly sorting can be slow + std::make_heap(e_ngb_later_begin, e_ngb_later_end, cmp1); + heapified=true; + } + // go back to storing an iterator in e_ngb_later? + time = events.lower_bound(e_ngb_later_begin->first); + still_dominated = true; + while (e_ngb_later_begin != e_ngb_later_end && e_ngb_later_begin->first <= time->first) { + Vertex w = e_ngb_later_begin->second; + auto wit = neighbors[dominator].find(w); + if (wit == neighbors[dominator].end() || wit->second > e_ngb_later_begin->first) + still_dominated = false; + e_ngb.insert(w); + std::pop_heap(e_ngb_later_begin, e_ngb_later_end--, cmp1); + } + } while (still_dominated); // this doesn't seem to help that much... + } +end_move: + if(dead) { + neighbors[u].erase(v); + neighbors[v].erase(u); + eg.erase(ei); + } else if(start_time != time){ + time->second.push_back(*ei); + neighbors[u][v]=time->first; + neighbors[v][u]=time->first; + eg.erase(ei); + } + } + // check if event is dead (all edges have moved) + if (evi->second.empty()) { + events.erase(evi); + } + } + } + + std::vector output() { + std::vector> r; + for(auto& ev : events) + for(auto& e : ev.second) + r.emplace_back(e.first, e.second, ev.first); + return r; + } + +}; + /** \brief Implicitly constructs a flag complex from edges as an input, collapses edges while preserving the persistent * homology and returns the remaining edges as a range. * @@ -377,6 +598,30 @@ template auto flag_complex_collapse_edges(const Filtere return remaining_edges; } +// Would it help to label the points according to some spatial sorting? +template auto flag_complex_collapse_edges2(FilteredEdgeRange&& edges, Delay&&delay) { + auto first_edge_itr = std::begin(edges); + using Vertex = std::decay_t(*first_edge_itr))>; + using Filtration_value = std::decay_t(*first_edge_itr))>; + using Edge_collapser = Flag_complex_edge_collapser2; + if (first_edge_itr != std::end(edges)) { + std::vector edges2; + if(std::is_same>::value) + edges2 = std::move(edges); + else + edges2.insert(edges2.end(), edges.begin(), edges.end()); + // The sorting is not necessary, but inserting in the map is faster if done in order + std::sort(edges2.begin(), edges2.end(), [](auto const&a, auto const&b){return std::get<2>(a)>std::get<2>(b);}); + Edge_collapser edge_collapser; + edge_collapser.process_edges(edges2, std::forward(delay)); + return edge_collapser.output(); + } + return std::vector(); +} +template auto flag_complex_collapse_edges2(const FilteredEdgeRange& edges) { + return flag_complex_collapse_edges2(edges, [](auto const&d){return d;}); +} + } // namespace collapse } // namespace Gudhi -- cgit v1.2.3