19#ifndef PERSISTENCE_ON_RECTANGLE_H
20#define PERSISTENCE_ON_RECTANGLE_H
22#include <gudhi/Debug_utils.h>
23#ifdef GUDHI_DETAILED_TIMES
24 #include <gudhi/Clock.h>
27#include <boost/range/adaptor/reversed.hpp>
30 #include <tbb/parallel_sort.h>
42namespace Gudhi::cubical_complex {
48template <
class Filtration_value,
class Index = std::
size_t,
bool output_index = false>
49struct Persistence_on_rectangle {
54 T_with_index() =
default;
56 bool operator<(T_with_index
const& other)
const {
57 return std::tie(first, second) < std::tie(other.first, other.second);
59 Index out()
const {
return second; }
64 T_no_index() =
default;
66 bool operator<(T_no_index
const& other)
const {
return first < other.first; }
69 typedef std::conditional_t<output_index, T_with_index, T_no_index> T;
75 Index size_x, size_y, input_size;
83 std::unique_ptr<T[]> data_v_;
84 T& data_vertex(Index i){
return data_v_[i]; }
85 T data_vertex(Index i)
const {
return data_v_[i]; }
87 std::conditional_t<output_index, Index, Filtration_value> global_min;
92 std::unique_ptr<Index[]> ds_parent_v_;
93 std::vector<Index> ds_parent_s_;
94 Index& ds_parent_vertex(Index n) {
return ds_parent_v_[n]; }
95 Index& ds_parent_square(Index n) {
return ds_parent_s_[n]; }
97 template<
class Parent>
98 Index ds_find_set_(Index v, Parent&&ds_parent) {
106 Index ancestor = ds_parent(v);
107 while (ancestor != v)
110 ancestor = ds_parent(v);
113 while (ancestor != v)
115 ds_parent(old) = ancestor;
122 Index parent = ds_parent(v);
123 Index grandparent = ds_parent(parent);
124 while (parent != grandparent)
126 ds_parent(v) = grandparent;
128 parent = ds_parent(v);
129 grandparent = ds_parent(parent);
134 Index parent = ds_parent(v);
135 Index grandparent = ds_parent(parent);
136 while (parent != grandparent)
138 ds_parent(v) = grandparent;
140 parent = grandparent;
141 grandparent = ds_parent(parent);
147 while (v != (parent = ds_parent(v)))
152 Index ds_find_set_vertex(Index v) {
153 return ds_find_set_(v, [
this](Index i) -> Index& {
return ds_parent_vertex(i); });
155 Index ds_find_set_square(Index v) {
156 return ds_find_set_(v, [
this](Index i) -> Index& {
return ds_parent_square(i); });
163 Edge(T f, Index v1, Index v2) : f(f), v1(v1), v2(v2) {}
165 bool operator<(Edge
const& other)
const {
return filt() < other.filt(); }
167 void dualize_edge(Edge& e)
const {
168 Index new_v2 = e.v1 + (dy + 1);
172 std::vector<Edge> edges;
175 input_size = n_rows * n_cols;
178 std::clog <<
"Input\n";
179 for(Index i = 0; i < input_size; ++i) {
180 std::clog << i <<
'\t' << input(i) <<
'\n';
187 data_v_.reset(
new T[input_size - dy - 1]);
188 ds_parent_v_.reset(
new Index[input_size - dy - 1]);
190 ds_parent_s_.resize(input_size);
192 edges.reserve(input_size / 2);
197 GUDHI_CHECK(a != b, std::logic_error(
"Bug in Gudhi: comparing a cell to itself"));
199 if (fb < fa)
return true;
200 if (fa < fb)
return false;
203 void set_parent_vertex(Index child, Index parent) {
204 GUDHI_CHECK(child != parent, std::logic_error(
"Bug in Gudhi: use mark_*_critical instead of set_parent"));
205 ds_parent_vertex(child) = parent;
207 void set_parent_square(Index child, Index parent) {
208 GUDHI_CHECK(child != parent, std::logic_error(
"Bug in Gudhi: use mark_*_critical instead of set_parent"));
209 ds_parent_square(child) = parent;
215 void fill_and_pair() {
218 auto mark_vertex_critical = [&](Index c) {
219 ds_parent_vertex(c) = c;
220 data_vertex(c) = T(f, i);
222 auto mark_square_critical = [&]() {
223 ds_parent_square(i) = i;
225 auto mark_edge_critical = [&](Index v1, Index v2) {
226 edges.emplace_back(T(f, i), v1, v2);
228 auto v_up_left = [&](){
return i - 1; };
229 auto v_up_right = [&](){
return i; };
230 auto v_down_left = [&](){
return i - dy - 1; };
231 auto v_down_right = [&](){
return i - dy; };
232 auto pair_square_up = [&](){ set_parent_square(i, i + dy); };
233 auto pair_square_down = [&](){ set_parent_square(i, i - dy); };
234 auto pair_square_left = [&](){ set_parent_square(i, i - 1); };
235 auto pair_square_right = [&](){ set_parent_square(i, i + 1); };
239 mark_vertex_critical(v_up_right());
240 i = size_x; f = input(i);
241 mark_vertex_critical(v_up_left());
242 i = dy * size_y; f = input(i);
243 mark_vertex_critical(v_down_right());
244 i = size_x + dy * size_y; f = input(i);
245 mark_vertex_critical(v_down_left());
248 for(Index x = 1; x < size_x; ++x) {
251 if (has_larger_input(i + dy, i, f)) {
252 auto up_left = [&](){
return has_larger_input(i - 1, i, f) && has_larger_input(i + dy - 1, i, f); };
253 auto up_right = [&](){
return has_larger_input(i + 1, i, f) && has_larger_input(i + dy + 1, i, f); };
255 set_parent_vertex(v_up_left(), v_up_right());
256 if (up_right()) mark_vertex_critical(v_up_right());
257 }
else if (up_right()) {
258 set_parent_vertex(v_up_right(), v_up_left());
260 mark_edge_critical(v_up_left(), v_up_right());
265 for(Index y = 1; y < size_y; ++y) {
270 if (has_larger_input(i + 1, i, f)) {
271 auto down_right = [&](){
return has_larger_input(i - dy, i, f) && has_larger_input(i + 1 - dy, i, f); };
272 auto up_right = [&](){
return has_larger_input(i + dy, i, f) && has_larger_input(i + 1 + dy, i, f); };
274 set_parent_vertex(v_down_right(), v_up_right());
275 if (up_right()) mark_vertex_critical(v_up_right());
276 }
else if (up_right()) {
277 set_parent_vertex(v_up_right(), v_down_right());
279 mark_edge_critical(v_down_right(), v_up_right());
284 for(Index x = 1; x < size_x; ++x) {
288 auto left = [&]() {
return has_larger_input(i - 1, i, f); };
289 auto right = [&]() {
return has_larger_input(i + 1, i, f); };
290 auto down = [&]() {
return has_larger_input(i - dy, i, f); };
291 auto up = [&]() {
return has_larger_input(i + dy, i, f); };
292 auto down_left = [&]() {
return has_larger_input(i - dy - 1, i, f); };
293 auto up_left = [&]() {
return has_larger_input(i + dy - 1, i, f); };
294 auto down_right = [&]() {
return has_larger_input(i - dy + 1, i, f); };
295 auto up_right = [&]() {
return has_larger_input(i + dy + 1, i, f); };
299 set_parent_vertex(v_up_left(), v_up_right());
302 set_parent_vertex(v_down_left(), v_up_left());
305 set_parent_vertex(v_down_right(), v_down_left());
308 mark_vertex_critical(v_up_right());
313 set_parent_vertex(v_up_right(), v_down_right());
315 mark_edge_critical(v_down_right(), v_up_right());
325 set_parent_vertex(v_down_right(), v_down_left());
327 mark_edge_critical(v_down_left(), v_down_right());
330 set_parent_vertex(v_up_right(), v_down_right());
332 mark_edge_critical(v_down_right(), v_up_right());
335 mark_edge_critical(v_down_left(), v_down_right());
342 set_parent_vertex(v_up_right(), v_down_right());
344 mark_edge_critical(v_down_right(), v_up_right());
352 set_parent_vertex(v_down_left(), v_up_left());
354 mark_edge_critical(v_down_left(), v_up_left());
358 set_parent_vertex(v_down_right(), v_down_left());
360 mark_edge_critical(v_down_left(), v_down_right());
363 set_parent_vertex(v_up_right(), v_down_right());
365 mark_edge_critical(v_down_right(), v_up_right());
368 mark_edge_critical(v_down_left(), v_down_right());
371 mark_edge_critical(v_down_left(), v_up_left());
374 set_parent_vertex(v_up_right(), v_down_right());
376 mark_edge_critical(v_down_right(), v_up_right());
386 set_parent_vertex(v_down_right(), v_down_left());
388 mark_edge_critical(v_down_left(), v_down_right());
391 set_parent_vertex(v_up_right(), v_down_right());
393 mark_edge_critical(v_down_right(), v_up_right());
396 mark_edge_critical(v_down_left(), v_down_right());
401 set_parent_vertex(v_up_right(), v_down_right());
403 mark_edge_critical(v_down_right(), v_up_right());
412 set_parent_vertex(v_down_left(), v_up_left());
415 set_parent_vertex(v_down_right(), v_down_left());
417 mark_edge_critical(v_down_left(), v_down_right());
427 set_parent_vertex(v_down_right(), v_down_left());
429 mark_edge_critical(v_down_left(), v_down_right());
431 mark_edge_critical(v_down_right(), v_up_right());
433 mark_edge_critical(v_down_left(), v_down_right());
439 mark_edge_critical(v_down_right(), v_up_right());
447 set_parent_vertex(v_down_right(), v_up_right());
449 mark_edge_critical(v_down_right(), v_up_right());
456 mark_square_critical();
466 if (has_larger_input(i - 1, i, f)) {
467 auto down_left = [&](){
return has_larger_input(i - dy, i, f) && has_larger_input(i - 1 - dy, i, f); };
468 auto up_left = [&](){
return has_larger_input(i + dy, i, f) && has_larger_input(i - 1 + dy, i, f); };
470 set_parent_vertex(v_down_left(), v_up_left());
471 if (up_left()) mark_vertex_critical(v_up_left());
472 }
else if (up_left()) {
473 set_parent_vertex(v_up_left(), v_down_left());
475 mark_edge_critical(v_down_left(), v_up_left());
481 for(Index x = 1; x < size_x; ++x) {
484 if (has_larger_input(i - dy, i, f)) {
485 auto down_left = [&](){
return has_larger_input(i - 1, i, f) && has_larger_input(i - dy - 1, i, f); };
486 auto down_right = [&](){
return has_larger_input(i + 1, i, f) && has_larger_input(i - dy + 1, i, f); };
488 set_parent_vertex(v_down_left(), v_down_right());
489 if (down_right()) mark_vertex_critical(v_down_right());
490 }
else if (down_right()) {
491 set_parent_vertex(v_down_right(), v_down_left());
493 mark_edge_critical(v_down_left(), v_down_right());
504 tbb::parallel_sort(edges.begin(), edges.end());
506 std::sort(edges.begin(), edges.end());
509 std::clog <<
"edges\n";
510 for(
auto&e : edges){ std::clog << e.v1 <<
'\t' << e.v2 <<
'\t' << e.filt() <<
'\n'; }
515 void primal(Out&&out){
516 auto it = std::remove_if(edges.begin(), edges.end(), [&](Edge& e) {
518 Index a = ds_find_set_vertex(e.v1);
519 Index b = ds_find_set_vertex(e.v2);
520 if (a == b) return false;
521 if (data_vertex(b) < data_vertex(a)) std::swap(a, b);
522 ds_parent_vertex(b) = a;
523 out(data_vertex(b).out(), e.f.out());
526 edges.erase(it, edges.end());
527 global_min = data_vertex(ds_find_set_vertex(0)).out();
534 for (
auto e : boost::adaptors::reverse(edges)) {
536 Index a = ds_find_set_square(e.v1);
537 Index b = ds_find_set_square(e.v2);
538 GUDHI_CHECK(a != b, std::logic_error(
"Bug in Gudhi"));
540 if (b == 0 || (a != 0 && input(a) < input(b))) std::swap(a, b);
541 ds_parent_square(b) = a;
542 if constexpr (output_index)
545 out(e.f.out(), input(b));
584template <
bool output_index = false,
typename Filtration_value,
typename Index,
typename Out0,
typename Out1>
585auto persistence_on_rectangle_from_top_cells(
Filtration_value const* input, Index n_rows, Index n_cols,
586 Out0&&out0, Out1&&out1){
587#ifdef GUDHI_DETAILED_TIMES
590 GUDHI_CHECK(n_rows >= 2 && n_cols >= 2,
591 std::domain_error(
"The complex must truly be 2d, i.e. at least 2 rows and 2 columns"));
592 Persistence_on_rectangle<Filtration_value, Index, output_index> X;
593 X.init(input, n_rows, n_cols);
594#ifdef GUDHI_DETAILED_TIMES
595 std::clog <<
"init: " << clock; clock.begin();
598#ifdef GUDHI_DETAILED_TIMES
599 std::clog <<
"fill and pair: " << clock; clock.begin();
602#ifdef GUDHI_DETAILED_TIMES
603 std::clog <<
"sort: " << clock; clock.begin();
606#ifdef GUDHI_DETAILED_TIMES
607 std::clog <<
"primal pass: " << clock; clock.begin();
610#ifdef GUDHI_DETAILED_TIMES
611 std::clog <<
"dual pass: " << clock;
Value type for a filtration function on a cell complex.
Definition FiltrationValue.h:20