diff options
Diffstat (limited to 'ripser.cpp')
| -rw-r--r-- | ripser.cpp | 89 |
1 files changed, 80 insertions, 9 deletions
@@ -36,12 +36,12 @@ */ -//#define USE_COEFFICIENTS +#define USE_COEFFICIENTS -//#define INDICATE_PROGRESS +#define INDICATE_PROGRESS #define PRINT_PERSISTENCE_PAIRS -//#define USE_ROBINHOOD_HASHMAP +#define USE_ROBINHOOD_HASHMAP #include <algorithm> #include <cassert> @@ -281,6 +281,14 @@ struct sparse_distance_matrix { } } + value_t operator()(const index_t i, const index_t j) const { + auto neighbor = + std::lower_bound(neighbors[i].begin(), neighbors[i].end(), index_diameter_t{j, 0}); + return (neighbor != neighbors[i].end() && get_index(*neighbor) == j) + ? get_diameter(*neighbor) + : std::numeric_limits<value_t>::infinity(); + } + size_t size() const { return neighbors.size(); } }; @@ -385,6 +393,7 @@ template <typename DistanceMatrix> class ripser { const binomial_coeff_table binomial_coeff; const std::vector<coefficient_t> multiplicative_inverse; mutable std::vector<diameter_entry_t> cofacet_entries; + mutable std::vector<index_t> vertices; struct entry_hash { std::size_t operator()(const entry_t& e) const { return hash<index_t>()(::get_index(e)); } @@ -426,6 +435,19 @@ public: return out; } + value_t compute_diameter(const index_t index, index_t dim) const { + value_t diam = -std::numeric_limits<value_t>::infinity(); + + vertices.clear(); + get_simplex_vertices(index, dim, dist.size(), std::back_inserter(vertices)); + + for (index_t i = 0; i <= dim; ++i) + for (index_t j = 0; j < i; ++j) { + diam = std::max(diam, dist(vertices[i], vertices[j])); + } + return diam; + } + class simplex_coboundary_enumerator; class simplex_boundary_enumerator { @@ -465,6 +487,41 @@ public: } }; + + diameter_entry_t get_apparent_facet(const diameter_entry_t simplex, const index_t dim) { + simplex_boundary_enumerator facets(simplex, dim, *this); + while (facets.has_next()) { + diameter_entry_t facet = facets.next(); + if (get_diameter(facet) == get_diameter(simplex)) { + simplex_coboundary_enumerator cofacets(facet, dim - 1, *this); + while (cofacets.has_next()) { + auto cofacet = cofacets.next(); + if (get_diameter(cofacet) == get_diameter(simplex)) + return (get_index(cofacet) == get_index(simplex)) ? facet + : std::make_pair(0, -1); + } + } + } + return std::make_pair(0, -1); + } + + diameter_entry_t get_apparent_cofacet(const diameter_entry_t simplex, const index_t dim) { + simplex_coboundary_enumerator cofacets(simplex, dim, *this); + while (cofacets.has_next()) { + diameter_entry_t cofacet = cofacets.next(); + if (get_diameter(cofacet) == get_diameter(simplex)) { + simplex_boundary_enumerator facets(cofacet, dim + 1, *this); + while (facets.has_next()) { + auto facet = facets.next(); + if (get_diameter(facet) == get_diameter(simplex)) + return (get_index(facet) == get_index(simplex)) ? cofacet + : std::make_pair(0, -1); + } + } + } + return std::make_pair(0, -1); + } + void assemble_columns_to_reduce(std::vector<diameter_index_t>& simplices, std::vector<diameter_index_t>& columns_to_reduce, entry_hash_map& pivot_column_index, index_t dim) { @@ -495,7 +552,9 @@ public: next_simplices.push_back({get_diameter(cofacet), get_index(cofacet)}); - if (pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()) + if ((pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()) && + (get_index(get_apparent_cofacet(cofacet, dim + 1)) == -1) && + (get_index(get_apparent_facet(cofacet, dim + 1)) == -1)) columns_to_reduce.push_back({get_diameter(cofacet), get_index(cofacet)}); } } @@ -537,8 +596,11 @@ public: std::cout << " [0," << get_diameter(e) << ")" << std::endl; #endif dset.link(u, v); - } else - columns_to_reduce.push_back(e); + } else { + if ((get_index(get_apparent_cofacet(e, 1)) == -1) && + (get_index(get_apparent_facet(e, 1)) == -1)) + columns_to_reduce.push_back(e); + } } std::reverse(columns_to_reduce.begin(), columns_to_reduce.end()); @@ -583,7 +645,7 @@ 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) { - bool check_for_emergent_pair = true; + bool check_for_emergent_pair = false; cofacet_entries.clear(); simplex_coboundary_enumerator cofacets(simplex, dim, *this); while (cofacets.has_next()) { @@ -651,8 +713,8 @@ public: greater_diameter_or_smaller_index<diameter_entry_t>> working_reduction_column, working_coboundary; - diameter_entry_t pivot = init_coboundary_and_get_pivot( - column_to_reduce, working_coboundary, dim, pivot_column_index); + diameter_entry_t e, pivot = init_coboundary_and_get_pivot( + column_to_reduce, working_coboundary, dim, pivot_column_index); while (true) { #ifdef INDICATE_PROGRESS @@ -677,6 +739,15 @@ public: factor, dim, working_reduction_column, working_coboundary); pivot = get_pivot(working_coboundary); + } else if (get_index(e = get_apparent_facet(pivot, dim + 1)) != -1) { + + set_coefficient(e, modulus - get_coefficient(e)); + + add_simplex_coboundary(e, dim, working_reduction_column, + working_coboundary); + + pivot = get_pivot(working_coboundary); + } else { #ifdef PRINT_PERSISTENCE_PAIRS value_t death = get_diameter(pivot); |