diff options
| -rw-r--r-- | Makefile | 4 | ||||
| -rw-r--r-- | ripser.cpp | 204 |
2 files changed, 161 insertions, 47 deletions
@@ -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 @@ -184,6 +184,9 @@ struct diameter_entry_t : std::pair<value_t, entry_t> { diameter_entry_t(const diameter_index_t& _diameter_index, coefficient_t _coefficient) : diameter_entry_t(get_diameter(_diameter_index), make_entry(get_index(_diameter_index), _coefficient)) {} + diameter_entry_t(const diameter_index_t& _diameter_index) + : diameter_entry_t(get_diameter(_diameter_index), + make_entry(get_index(_diameter_index), 0)) {} diameter_entry_t(const index_t& _index) : diameter_entry_t(0, _index, 0) {} }; @@ -198,13 +201,17 @@ void set_coefficient(diameter_entry_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); } -template <typename Entry> struct greater_diameter_or_smaller_index { +template <typename Entry> 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 <typename Entry> 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 <compressed_matrix_layout Layout> struct compressed_distance_matrix { @@ -265,9 +272,6 @@ struct sparse_distance_matrix { index_t num_edges; - mutable std::vector<std::vector<index_diameter_t>::const_reverse_iterator> neighbor_it; - mutable std::vector<std::vector<index_diameter_t>::const_reverse_iterator> neighbor_end; - sparse_distance_matrix(std::vector<std::vector<index_diameter_t>>&& _neighbors, index_t _num_edges) : neighbors(std::move(_neighbors)), num_edges(_num_edges) {} @@ -284,6 +288,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(); } }; @@ -388,6 +400,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)); } @@ -429,15 +442,28 @@ public: return out; } + value_t compute_diameter(const index_t index, const index_t dim) const { + value_t diam = -std::numeric_limits<value_t>::infinity(); + + vertices.resize(dim + 1); + get_simplex_vertices(index, dim, dist.size(), vertices.rbegin()); + + 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 { private: index_t idx_below, idx_above, j, k; - const diameter_entry_t simplex; + diameter_entry_t simplex; + index_t dim; const coefficient_t modulus; const binomial_coeff_table& binomial_coeff; - const index_t dim; const ripser& parent; public: @@ -445,7 +471,19 @@ public: 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) {} + 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; + dim = _dim; + } bool has_next() { return (k >= 0); } @@ -468,6 +506,47 @@ public: } }; + diameter_entry_t get_zero_pivot_facet(const diameter_entry_t simplex, const index_t dim) { + static simplex_boundary_enumerator facets(0, *this); + facets.set_simplex(simplex, dim); + while (facets.has_next()) { + diameter_entry_t facet = facets.next(); + if (get_diameter(facet) == get_diameter(simplex)) return facet; + } + return diameter_entry_t(-1); + } + + diameter_entry_t get_zero_pivot_cofacet(const diameter_entry_t simplex, const index_t dim) { + static simplex_coboundary_enumerator cofacets(*this); + cofacets.set_simplex(simplex, dim); + while (cofacets.has_next()) { + diameter_entry_t cofacet = cofacets.next(); + if (get_diameter(cofacet) == get_diameter(simplex)) return cofacet; + } + return diameter_entry_t(-1); + } + + diameter_entry_t get_zero_apparent_facet(const diameter_entry_t simplex, const index_t dim) { + diameter_entry_t facet = get_zero_pivot_facet(simplex, dim); + return ((get_index(facet) != -1) && + (get_index(get_zero_pivot_cofacet(facet, dim - 1)) == get_index(simplex))) + ? facet + : diameter_entry_t(-1); + } + + diameter_entry_t get_zero_apparent_cofacet(const diameter_entry_t simplex, const index_t dim) { + diameter_entry_t cofacet = get_zero_pivot_cofacet(simplex, dim); + return ((get_index(cofacet) != -1) && + (get_index(get_zero_pivot_facet(cofacet, dim + 1)) == get_index(simplex))) + ? cofacet + : diameter_entry_t(-1); + } + + bool is_in_zero_apparent_pair(const diameter_entry_t simplex, const index_t dim) { + return (get_index(get_zero_apparent_cofacet(simplex, dim)) != -1) || + (get_index(get_zero_apparent_facet(simplex, dim)) != -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) { @@ -477,12 +556,13 @@ public: std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step; #endif - --dim; columns_to_reduce.clear(); std::vector<diameter_index_t> next_simplices; + simplex_coboundary_enumerator cofacets(*this); + for (diameter_index_t& simplex : simplices) { - simplex_coboundary_enumerator cofacets(diameter_entry_t(simplex, 1), dim, *this); + cofacets.set_simplex(diameter_entry_t(simplex, 1), dim - 1); while (cofacets.has_next(false)) { #ifdef INDICATE_PROGRESS @@ -495,10 +575,9 @@ public: #endif auto cofacet = cofacets.next(); if (get_diameter(cofacet) <= threshold) { - next_simplices.push_back({get_diameter(cofacet), get_index(cofacet)}); - - if (pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()) + if (!is_in_zero_apparent_pair(cofacet, dim) && + (pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end())) columns_to_reduce.push_back({get_diameter(cofacet), get_index(cofacet)}); } } @@ -512,7 +591,7 @@ public: #endif std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), - greater_diameter_or_smaller_index<diameter_index_t>()); + greater_diameter_or_smaller_index<diameter_index_t>); #ifdef INDICATE_PROGRESS std::cerr << clear_line << std::flush; #endif @@ -528,7 +607,7 @@ public: edges = get_edges(); std::sort(edges.rbegin(), edges.rend(), - greater_diameter_or_smaller_index<diameter_index_t>()); + greater_diameter_or_smaller_index<diameter_index_t>); std::vector<index_t> vertices_of_edge(2); for (auto e : edges) { get_simplex_vertices(get_index(e), 1, n, vertices_of_edge.rbegin()); @@ -540,7 +619,7 @@ public: std::cout << " [0," << get_diameter(e) << ")" << std::endl; #endif dset.link(u, v); - } else + } else if (get_index(get_zero_apparent_cofacet(e, 1)) == -1) columns_to_reduce.push_back(e); } std::reverse(columns_to_reduce.begin(), columns_to_reduce.end()); @@ -586,15 +665,17 @@ 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) { + static simplex_coboundary_enumerator cofacets(*this); bool check_for_emergent_pair = true; cofacet_entries.clear(); - simplex_coboundary_enumerator cofacets(simplex, dim, *this); + cofacets.set_simplex(simplex, dim); while (cofacets.has_next()) { diameter_entry_t cofacet = cofacets.next(); if (get_diameter(cofacet) <= threshold) { 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()) + if ((pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()) && + (get_index(get_zero_apparent_facet(cofacet, dim + 1)) == -1)) return cofacet; check_for_emergent_pair = false; } @@ -607,8 +688,9 @@ public: template <typename Column> void add_simplex_coboundary(const diameter_entry_t simplex, const index_t& dim, Column& working_reduction_column, Column& working_coboundary) { + static simplex_coboundary_enumerator cofacets(*this); working_reduction_column.push(simplex); - simplex_coboundary_enumerator cofacets(simplex, dim, *this); + cofacets.set_simplex(simplex, dim); while (cofacets.has_next()) { diameter_entry_t cofacet = cofacets.next(); if (get_diameter(cofacet) <= threshold) working_coboundary.push(cofacet); @@ -638,7 +720,7 @@ public: #endif compressed_sparse_matrix<diameter_entry_t> reduction_matrix; - + #ifdef INDICATE_PROGRESS std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step; #endif @@ -651,11 +733,11 @@ public: reduction_matrix.append_column(); std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, - greater_diameter_or_smaller_index<diameter_entry_t>> + greater_diameter_or_smaller_index_comp<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 @@ -680,6 +762,12 @@ public: factor, dim, working_reduction_column, working_coboundary); pivot = get_pivot(working_coboundary); + } else if (get_index(e = get_zero_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); @@ -739,17 +827,31 @@ public: template <> class ripser<compressed_lower_distance_matrix>::simplex_coboundary_enumerator { index_t idx_below, idx_above, j, k; std::vector<index_t> vertices; - const 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; 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) { + const ripser& _parent) + : 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()); + } + + simplex_coboundary_enumerator(const ripser& _parent) : modulus(_parent.modulus), dist(_parent.dist), + binomial_coeff(_parent.binomial_coeff), parent(_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; + simplex = _simplex; + vertices.resize(_dim + 1); parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin()); } @@ -777,28 +879,41 @@ public: template <> class ripser<sparse_distance_matrix>::simplex_coboundary_enumerator { index_t idx_below, idx_above, k; std::vector<index_t> vertices; - const 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<std::vector<index_diameter_t>::const_reverse_iterator>& neighbor_it; - std::vector<std::vector<index_diameter_t>::const_reverse_iterator>& neighbor_end; + std::vector<std::vector<index_diameter_t>::const_reverse_iterator> neighbor_it; + std::vector<std::vector<index_diameter_t>::const_reverse_iterator> neighbor_end; index_diameter_t neighbor; + 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(dist.neighbor_it), - neighbor_end(dist.neighbor_end) { - neighbor_it.clear(); - neighbor_end.clear(); - _parent.get_simplex_vertices(idx_below, _dim, _parent.n, vertices.rbegin()); - - for (auto v : vertices) { - neighbor_it.push_back(dist.neighbors[v].rbegin()); - neighbor_end.push_back(dist.neighbors[v].rend()); + : modulus(_parent.modulus), dist(_parent.dist), + binomial_coeff(_parent.binomial_coeff), parent(_parent) { + if (get_index(_simplex) != -1) set_simplex(_simplex, _dim); + } + + simplex_coboundary_enumerator(const ripser& _parent) + : modulus(_parent.modulus), dist(_parent.dist), + binomial_coeff(_parent.binomial_coeff), parent(_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; + vertices.resize(_dim + 1); + parent.get_simplex_vertices(idx_below, _dim, parent.n, vertices.rbegin()); + + neighbor_it.resize(_dim + 1); + neighbor_end.resize(_dim + 1); + 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(); } } @@ -1145,8 +1260,7 @@ int main(int argc, char** argv) { for (auto d : dist.distances) { min = std::min(min, d); max = std::max(max, d); - max_finite = - d != std::numeric_limits<value_t>::infinity() ? std::max(max, d) : max_finite; + if (d != std::numeric_limits<value_t>::infinity()) max_finite = std::max(max_finite, d); if (d <= threshold) ++num_edges; } std::cout << "value range: [" << min << "," << max_finite << "]" << std::endl; |