From fd53851baabfd89b23c620826a8b6265418951a2 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 27 Jan 2021 20:58:24 +0100 Subject: restructured apparent pairs --- Makefile | 4 +- ripser.cpp | 286 +++++++++++++++++++++++++++++++------------------------------ 2 files changed, 147 insertions(+), 143 deletions(-) diff --git a/Makefile b/Makefile index ab410bd..921e438 100644 --- a/Makefile +++ b/Makefile @@ -5,10 +5,10 @@ all: ripser ripser-coeff ripser-debug ripser: ripser.cpp - c++ -std=c++11 -Wall ripser.cpp -o ripser -Ofast -D NDEBUG + c++ -std=c++11 -Wall ripser.cpp -o ripser -O3 -D NDEBUG ripser-coeff: ripser.cpp - c++ -std=c++11 -Wall ripser.cpp -o ripser-coeff -Ofast -D NDEBUG -D USE_COEFFICIENTS + c++ -std=c++11 -Wall ripser.cpp -o ripser-coeff -O3 -D NDEBUG -D USE_COEFFICIENTS ripser-debug: ripser.cpp c++ -std=c++11 -Wall ripser.cpp -o ripser-debug -g diff --git a/ripser.cpp b/ripser.cpp index 1ffa443..6f0fd0a 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -201,13 +201,17 @@ void set_coefficient(diameter_entry_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); } -template struct greater_diameter_or_smaller_index { +template struct greater_diameter_or_smaller_index_comp { bool operator()(const Entry& a, const Entry& b) { - return (get_diameter(a) > get_diameter(b)) || - ((get_diameter(a) == get_diameter(b)) && (get_index(a) < get_index(b))); + return greater_diameter_or_smaller_index(a, b); } }; +template bool greater_diameter_or_smaller_index(const Entry& a, const Entry& b) { + return (get_diameter(a) > get_diameter(b)) || + ((get_diameter(a) == get_diameter(b)) && (get_index(a) < get_index(b))); +} + enum compressed_matrix_layout { LOWER_TRIANGULAR, UPPER_TRIANGULAR }; template struct compressed_distance_matrix { @@ -456,32 +460,31 @@ public: class simplex_boundary_enumerator { private: index_t idx_below, idx_above, j, k; - diameter_entry_t simplex; + diameter_entry_t simplex; const coefficient_t modulus; const binomial_coeff_table& binomial_coeff; const index_t dim; const ripser& parent; public: - simplex_boundary_enumerator(const diameter_entry_t _simplex, const index_t _dim, - const ripser& _parent) - : idx_below(get_index(_simplex)), idx_above(0), j(_parent.n - 1), k(_dim), - simplex(_simplex), modulus(_parent.modulus), binomial_coeff(_parent.binomial_coeff), - dim(_dim), parent(_parent) {} - - simplex_boundary_enumerator(const index_t _dim, const ripser& _parent) - : simplex_boundary_enumerator(-1, _dim, _parent) {} - - void set_simplex(const diameter_entry_t _simplex, const index_t _dim) { - idx_below = get_index(_simplex); - idx_above = 0; - j = parent.n - 1; - k = _dim; - simplex = _simplex; - } - - - bool has_next() { return (k >= 0); } + simplex_boundary_enumerator(const diameter_entry_t _simplex, const index_t _dim, + const ripser& _parent) + : idx_below(get_index(_simplex)), idx_above(0), j(_parent.n - 1), k(_dim), + simplex(_simplex), modulus(_parent.modulus), binomial_coeff(_parent.binomial_coeff), + dim(_dim), parent(_parent) {} + + simplex_boundary_enumerator(const index_t _dim, const ripser& _parent) + : simplex_boundary_enumerator(-1, _dim, _parent) {} + + void set_simplex(const diameter_entry_t _simplex, const index_t _dim) { + idx_below = get_index(_simplex); + idx_above = 0; + j = parent.n - 1; + k = _dim; + simplex = _simplex; + } + + bool has_next() { return (k >= 0); } diameter_entry_t next() { j = parent.get_max_vertex(idx_below, k + 1, j); @@ -502,57 +505,56 @@ public: } }; - diameter_entry_t get_apparent_facet(const diameter_entry_t simplex, const index_t dim, - simplex_boundary_enumerator& facets, - simplex_coboundary_enumerator& cofacets) { - facets.set_simplex(simplex, dim); + diameter_entry_t get_zero_pivot_facet(const diameter_entry_t simplex, const index_t dim, + simplex_boundary_enumerator& facets) { + facets.set_simplex(simplex, dim); while (facets.has_next()) { diameter_entry_t facet = facets.next(); - if (get_diameter(facet) == get_diameter(simplex)) { - cofacets.set_simplex(facet, dim - 1); - while (cofacets.has_next()) { - diameter_entry_t cofacet = cofacets.next(); - if (get_diameter(cofacet) == get_diameter(simplex)) - return (get_index(cofacet) == get_index(simplex)) ? facet - : diameter_entry_t(-1); - } - } + if (get_diameter(facet) == get_diameter(simplex)) return facet; } return diameter_entry_t(-1); } - diameter_entry_t get_apparent_facet(const diameter_entry_t simplex, const index_t dim) { - simplex_boundary_enumerator facets(dim, *this); - simplex_coboundary_enumerator cofacets(dim - 1, *this); - return get_apparent_facet(simplex, dim, facets, cofacets); - } - - diameter_entry_t get_apparent_cofacet(const diameter_entry_t simplex, const index_t dim, - simplex_coboundary_enumerator& cofacets, - simplex_boundary_enumerator& facets) { - cofacets.set_simplex(simplex, dim); + diameter_entry_t get_zero_pivot_cofacet(const diameter_entry_t simplex, const index_t dim, + simplex_coboundary_enumerator& cofacets) { + cofacets.set_simplex(simplex, dim); while (cofacets.has_next()) { - diameter_entry_t cofacet = cofacets.next(); - if (get_diameter(cofacet) == get_diameter(simplex)) { - - facets.set_simplex(cofacet, dim + 1); - while (facets.has_next()) { - auto facet = facets.next(); - if (get_diameter(facet) == get_diameter(simplex)) - return (get_index(facet) == get_index(simplex)) ? cofacet - : diameter_entry_t(-1); - } - } + diameter_entry_t cofacet = cofacets.next(); + if (get_diameter(cofacet) == get_diameter(simplex)) return cofacet; } return diameter_entry_t(-1); } - diameter_entry_t get_apparent_cofacet(const diameter_entry_t simplex, const index_t dim) { - simplex_coboundary_enumerator cofacets(dim, *this); - simplex_boundary_enumerator facets(dim + 1, *this); - return get_apparent_cofacet(simplex, dim, cofacets, facets); - } - + diameter_entry_t get_zero_apparent_facet(const diameter_entry_t simplex, const index_t dim, + simplex_boundary_enumerator& facets, + simplex_coboundary_enumerator& cofacets) { + diameter_entry_t facet = get_zero_pivot_facet(simplex, dim, facets); + return ((get_index(facet) != -1) && + (get_index(get_zero_pivot_cofacet(facet, dim - 1, cofacets)) == get_index(simplex))) + ? facet + : diameter_entry_t(-1); + } + + diameter_entry_t get_zero_apparent_cofacet(const diameter_entry_t simplex, const index_t dim, + simplex_coboundary_enumerator& cofacets, + simplex_boundary_enumerator& facets) { + diameter_entry_t cofacet = get_zero_pivot_cofacet(simplex, dim, cofacets); + return ((get_index(cofacet) != -1) && + (get_index(get_zero_pivot_facet(cofacet, dim + 1, facets)) == get_index(simplex))) + ? cofacet + : diameter_entry_t(-1); + } + + bool is_in_zero_apparent_pair(const diameter_entry_t simplex, const index_t dim, + simplex_coboundary_enumerator& cofacets_dim, + simplex_coboundary_enumerator& cofacets_dim_1, + simplex_boundary_enumerator& facets_dim, + simplex_boundary_enumerator& facets_dim_1) { + return (get_index(get_zero_apparent_cofacet(simplex, dim, cofacets_dim, facets_dim_1)) != + -1) || + (get_index(get_zero_apparent_facet(simplex, dim, facets_dim, cofacets_dim_1)) != -1); + } + void assemble_columns_to_reduce(std::vector& simplices, std::vector& columns_to_reduce, entry_hash_map& pivot_column_index, index_t dim) { @@ -564,9 +566,10 @@ public: columns_to_reduce.clear(); std::vector next_simplices; - - simplex_coboundary_enumerator cofacets(dim - 1, *this), cofacets_dim(dim, *this), cofacets_dim_1(dim - 1, *this); - simplex_boundary_enumerator facets_dim(dim, *this), facets_dim_1(dim + 1, *this); + + simplex_coboundary_enumerator cofacets(dim - 1, *this), cofacets_dim(dim, *this), + cofacets_dim_1(dim - 1, *this); + simplex_boundary_enumerator facets_dim(dim, *this), facets_dim_1(dim + 1, *this); for (diameter_index_t& simplex : simplices) { cofacets.set_simplex(diameter_entry_t(simplex, 1), dim - 1); @@ -582,12 +585,10 @@ public: #endif auto cofacet = cofacets.next(); if (get_diameter(cofacet) <= threshold) { - next_simplices.push_back({get_diameter(cofacet), get_index(cofacet)}); - - if ((get_index(get_apparent_cofacet(cofacet, dim, cofacets_dim, facets_dim_1)) == -1) && - (pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()) && - (get_index(get_apparent_facet(cofacet, dim, facets_dim, cofacets_dim_1)) == -1)) + if (!is_in_zero_apparent_pair(cofacet, dim, cofacets_dim, cofacets_dim_1, + facets_dim, facets_dim_1) && + (pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end())) columns_to_reduce.push_back({get_diameter(cofacet), get_index(cofacet)}); } } @@ -601,7 +602,7 @@ public: #endif std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), - greater_diameter_or_smaller_index()); + greater_diameter_or_smaller_index); #ifdef INDICATE_PROGRESS std::cerr << clear_line << std::flush; #endif @@ -614,13 +615,13 @@ public: #endif union_find dset(n); - + simplex_coboundary_enumerator cofacets(1, *this); simplex_boundary_enumerator facets(2, *this); edges = get_edges(); std::sort(edges.rbegin(), edges.rend(), - greater_diameter_or_smaller_index()); + greater_diameter_or_smaller_index); std::vector vertices_of_edge(2); for (auto e : edges) { get_simplex_vertices(get_index(e), 1, n, vertices_of_edge.rbegin()); @@ -632,7 +633,7 @@ public: std::cout << " [0," << get_diameter(e) << ")" << std::endl; #endif dset.link(u, v); - } else if (get_index(get_apparent_cofacet(e, 1, cofacets, facets)) == -1) + } else if (get_index(get_zero_apparent_cofacet(e, 1, cofacets, facets)) == -1) columns_to_reduce.push_back(e); } std::reverse(columns_to_reduce.begin(), columns_to_reduce.end()); @@ -678,10 +679,9 @@ public: diameter_entry_t init_coboundary_and_get_pivot(const diameter_entry_t simplex, Column& working_coboundary, const index_t& dim, entry_hash_map& pivot_column_index, - simplex_coboundary_enumerator& cofacets, - simplex_boundary_enumerator& facets, - simplex_coboundary_enumerator& other_cofacets - ) { + simplex_coboundary_enumerator& cofacets, + simplex_boundary_enumerator& facets, + simplex_coboundary_enumerator& other_cofacets) { bool check_for_emergent_pair = true; cofacet_entries.clear(); cofacets.set_simplex(simplex, dim); @@ -691,7 +691,8 @@ public: cofacet_entries.push_back(cofacet); if (check_for_emergent_pair && (get_diameter(simplex) == get_diameter(cofacet))) { if ((pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()) && - (get_index(get_apparent_facet(cofacet, dim + 1, facets, other_cofacets)) == -1)) + (get_index(get_zero_apparent_facet(cofacet, dim + 1, facets, + other_cofacets)) == -1)) return cofacet; check_for_emergent_pair = false; } @@ -704,9 +705,9 @@ public: template void add_simplex_coboundary(const diameter_entry_t simplex, const index_t& dim, Column& working_reduction_column, Column& working_coboundary, - simplex_coboundary_enumerator& cofacets) { + simplex_coboundary_enumerator& cofacets) { working_reduction_column.push(simplex); - cofacets.set_simplex(simplex, dim); + cofacets.set_simplex(simplex, dim); while (cofacets.has_next()) { diameter_entry_t cofacet = cofacets.next(); if (get_diameter(cofacet) <= threshold) working_coboundary.push(cofacet); @@ -718,14 +719,15 @@ public: const std::vector& columns_to_reduce, const size_t index_column_to_add, const coefficient_t factor, const size_t& dim, Column& working_reduction_column, - Column& working_coboundary, - simplex_coboundary_enumerator& cofacets) { + Column& working_coboundary, simplex_coboundary_enumerator& cofacets) { diameter_entry_t column_to_add(columns_to_reduce[index_column_to_add], factor); - add_simplex_coboundary(column_to_add, dim, working_reduction_column, working_coboundary, cofacets); + add_simplex_coboundary(column_to_add, dim, working_reduction_column, working_coboundary, + cofacets); for (diameter_entry_t simplex : reduction_matrix.subrange(index_column_to_add)) { set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); - add_simplex_coboundary(simplex, dim, working_reduction_column, working_coboundary, cofacets); + add_simplex_coboundary(simplex, dim, working_reduction_column, working_coboundary, + cofacets); } } @@ -737,10 +739,9 @@ public: #endif compressed_sparse_matrix reduction_matrix; - - simplex_coboundary_enumerator cofacets(dim, *this), other_cofacets(dim, *this); - simplex_boundary_enumerator facets(dim + 1, *this); + simplex_coboundary_enumerator cofacets(dim, *this), other_cofacets(dim, *this); + simplex_boundary_enumerator facets(dim + 1, *this); #ifdef INDICATE_PROGRESS std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step; @@ -754,11 +755,12 @@ public: reduction_matrix.append_column(); std::priority_queue, - greater_diameter_or_smaller_index> + greater_diameter_or_smaller_index_comp> working_reduction_column, working_coboundary; diameter_entry_t e, pivot = init_coboundary_and_get_pivot( - column_to_reduce, working_coboundary, dim, pivot_column_index, cofacets, facets, other_cofacets); + column_to_reduce, working_coboundary, dim, pivot_column_index, + cofacets, facets, other_cofacets); while (true) { #ifdef INDICATE_PROGRESS @@ -780,15 +782,17 @@ public: modulus; add_coboundary(reduction_matrix, columns_to_reduce, index_column_to_add, - factor, dim, working_reduction_column, working_coboundary, cofacets); + factor, dim, working_reduction_column, working_coboundary, + cofacets); pivot = get_pivot(working_coboundary); - } else if (get_index(e = get_apparent_facet(pivot, dim + 1, facets, other_cofacets)) != -1) { + } else if (get_index(e = get_zero_apparent_facet(pivot, dim + 1, facets, + other_cofacets)) != -1) { set_coefficient(e, modulus - get_coefficient(e)); - add_simplex_coboundary(e, dim, working_reduction_column, - working_coboundary, other_cofacets); + add_simplex_coboundary(e, dim, working_reduction_column, working_coboundary, + other_cofacets); pivot = get_pivot(working_coboundary); @@ -851,34 +855,34 @@ public: template <> class ripser::simplex_coboundary_enumerator { index_t idx_below, idx_above, j, k; std::vector vertices; - diameter_entry_t simplex; + diameter_entry_t simplex; const coefficient_t modulus; const compressed_lower_distance_matrix& dist; const binomial_coeff_table& binomial_coeff; - const ripser& parent; + const ripser& parent; public: simplex_coboundary_enumerator(const diameter_entry_t _simplex, const index_t _dim, const ripser& _parent) : idx_below(get_index(_simplex)), idx_above(0), j(_parent.n - 1), k(_dim + 1), vertices(_dim + 1), simplex(_simplex), modulus(_parent.modulus), dist(_parent.dist), - binomial_coeff(_parent.binomial_coeff), parent(_parent) { - if (get_index(_simplex) != -1) - parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin()); + binomial_coeff(_parent.binomial_coeff), parent(_parent) { + if (get_index(_simplex) != -1) + parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin()); + } + + simplex_coboundary_enumerator(const index_t _dim, const ripser& _parent) + : simplex_coboundary_enumerator(-1, _dim, _parent) {} + + void set_simplex(const diameter_entry_t _simplex, const index_t _dim) { + idx_below = get_index(_simplex); + idx_above = 0; + j = parent.n - 1; + k = _dim + 1; + vertices.resize(_dim + 1); + simplex = _simplex; + parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin()); } - - simplex_coboundary_enumerator(const index_t _dim, const ripser& _parent) - : simplex_coboundary_enumerator(-1, _dim, _parent) {} - - void set_simplex(const diameter_entry_t _simplex, const index_t _dim) { - idx_below = get_index(_simplex); - idx_above = 0; - j = parent.n - 1; - k = _dim + 1; - vertices.resize(_dim + 1); - simplex = _simplex; - parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin()); - } bool has_next(bool all_cofacets = true) { return (j >= k && (all_cofacets || binomial_coeff(j, k) > idx_below)); @@ -904,42 +908,42 @@ public: template <> class ripser::simplex_coboundary_enumerator { index_t idx_below, idx_above, k; std::vector vertices; - diameter_entry_t simplex; + diameter_entry_t simplex; const coefficient_t modulus; const sparse_distance_matrix& dist; const binomial_coeff_table& binomial_coeff; std::vector::const_reverse_iterator> neighbor_it; std::vector::const_reverse_iterator> neighbor_end; index_diameter_t neighbor; - const ripser& parent; + const ripser& parent; public: - simplex_coboundary_enumerator(const diameter_entry_t _simplex, const index_t _dim, - const ripser& _parent) - : idx_below(get_index(_simplex)), idx_above(0), k(_dim + 1), vertices(_dim + 1), - simplex(_simplex), modulus(_parent.modulus), dist(_parent.dist), - binomial_coeff(_parent.binomial_coeff), - neighbor_it(_dim + 1), neighbor_end(_dim + 1), parent(_parent) { - if (get_index(_simplex) != -1) set_simplex(_simplex, _dim); - } - - simplex_coboundary_enumerator(const index_t _dim, const ripser& _parent) - : simplex_coboundary_enumerator(-1, _dim, _parent) {} - - void set_simplex(const diameter_entry_t _simplex, const index_t _dim) { - idx_below = get_index(_simplex); - idx_above = 0; - k = _dim + 1; - simplex = _simplex; - - parent.get_simplex_vertices(idx_below, _dim, parent.n, vertices.rbegin()); - - for (index_t i = 0; i <= _dim; ++i) { - auto v = vertices[i]; - neighbor_it[i] = dist.neighbors[v].rbegin(); - neighbor_end[i] = dist.neighbors[v].rend(); - } - } + simplex_coboundary_enumerator(const diameter_entry_t _simplex, const index_t _dim, + const ripser& _parent) + : idx_below(get_index(_simplex)), idx_above(0), k(_dim + 1), vertices(_dim + 1), + simplex(_simplex), modulus(_parent.modulus), dist(_parent.dist), + binomial_coeff(_parent.binomial_coeff), neighbor_it(_dim + 1), neighbor_end(_dim + 1), + parent(_parent) { + if (get_index(_simplex) != -1) set_simplex(_simplex, _dim); + } + + simplex_coboundary_enumerator(const index_t _dim, const ripser& _parent) + : simplex_coboundary_enumerator(-1, _dim, _parent) {} + + void set_simplex(const diameter_entry_t _simplex, const index_t _dim) { + idx_below = get_index(_simplex); + idx_above = 0; + k = _dim + 1; + simplex = _simplex; + + parent.get_simplex_vertices(idx_below, _dim, parent.n, vertices.rbegin()); + + for (index_t i = 0; i <= _dim; ++i) { + auto v = vertices[i]; + neighbor_it[i] = dist.neighbors[v].rbegin(); + neighbor_end[i] = dist.neighbors[v].rend(); + } + } bool has_next(bool all_cofacets = true) { for (auto &it0 = neighbor_it[0], &end0 = neighbor_end[0]; it0 != end0; ++it0) { -- cgit v1.2.3