diff options
| -rw-r--r-- | ripser.cpp | 211 |
1 files changed, 171 insertions, 40 deletions
@@ -36,9 +36,9 @@ */ -//#define USE_COEFFICIENTS +#define USE_COEFFICIENTS -//#define INDICATE_PROGRESS +#define INDICATE_PROGRESS #define PRINT_PERSISTENCE_PAIRS //#define USE_GOOGLE_HASHMAP @@ -263,9 +263,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) {} @@ -281,6 +278,11 @@ struct sparse_distance_matrix { neighbors[i].push_back({j, mat(i, j)}); } } + + 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(); } }; @@ -386,6 +388,8 @@ 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 { @@ -429,7 +433,93 @@ 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 { + private: + index_t idx_below, idx_above, v, k; + const 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), v(_parent.n - 1), k(_dim + 1), + simplex(_simplex), modulus(_parent.modulus), + binomial_coeff(_parent.binomial_coeff), dim(_dim), parent(_parent) { + } + + bool has_next() { + return (k > 0); + } + + diameter_entry_t next() { + v = parent.get_max_vertex(idx_below, k, v); + + index_t face_index = idx_above - binomial_coeff(v, k) + idx_below; + + value_t face_diameter = parent.compute_diameter(face_index, dim - 1); + + coefficient_t face_coefficient = + (k & 1 ? 1 : -1 + modulus) * get_coefficient(simplex) % modulus; + + idx_below -= binomial_coeff(v, k); + idx_above += binomial_coeff(v, k - 1); + + --k; + + return diameter_entry_t(face_diameter, face_index, face_coefficient); + } + }; + + 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, @@ -461,7 +551,10 @@ 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)}); } } @@ -503,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()); @@ -549,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()) { @@ -618,7 +714,7 @@ public: diameter_entry_t pivot = init_coboundary_and_get_pivot( column_to_reduce, working_coboundary, dim, pivot_column_index); - + while (true) { #ifdef INDICATE_PROGRESS if (std::chrono::steady_clock::now() > next) { @@ -633,6 +729,7 @@ public: if (pair != pivot_column_index.end()) { entry_t other_pivot = pair->first; index_t index_column_to_add = pair->second; + coefficient_t factor = modulus - get_coefficient(pivot) * multiplicative_inverse[get_coefficient(other_pivot)] % @@ -643,24 +740,37 @@ public: pivot = get_pivot(working_coboundary); } else { + auto check = get_apparent_facet(pivot, dim + 1); + + if (get_index(check) != -1) { + + set_coefficient(check, modulus - get_coefficient(check)); + + add_simplex_coboundary(check, dim, working_reduction_column, + working_coboundary); + + pivot = get_pivot(working_coboundary); + + } else { #ifdef PRINT_PERSISTENCE_PAIRS - value_t death = get_diameter(pivot); - if (death > diameter * ratio) { + value_t death = get_diameter(pivot); + if (death > diameter * ratio) { #ifdef INDICATE_PROGRESS - std::cerr << clear_line << std::flush; + std::cerr << clear_line << std::flush; #endif - std::cout << " [" << diameter << "," << death << ")" << std::endl; - } + std::cout << " [" << diameter << "," << death << ")" << std::endl; + } #endif - pivot_column_index.insert({get_entry(pivot), index_column_to_reduce}); - - while (true) { - diameter_entry_t e = pop_pivot(working_reduction_column); - if (get_index(e) == -1) break; - assert(get_coefficient(e) > 0); - reduction_matrix.push_back(e); + pivot_column_index.insert({get_entry(pivot), index_column_to_reduce}); + + while (true) { + diameter_entry_t e = pop_pivot(working_reduction_column); + if (get_index(e) == -1) break; + assert(get_coefficient(e) > 0); + reduction_matrix.push_back(e); + } + break; } - break; } } else { #ifdef PRINT_PERSISTENCE_PAIRS @@ -701,17 +811,28 @@ public: template <> class ripser<compressed_lower_distance_matrix>::simplex_coboundary_enumerator { index_t idx_below, idx_above, v, 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), v(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) + : idx_below(get_index(_simplex)), idx_above(0), v(_parent.n - 1), k(_dim + 1), + vertices(_dim + 1), simplex(_simplex), modulus(_parent.modulus), dist(_parent.dist), + binomial_coeff(_parent.binomial_coeff), parent(_parent) { + parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin()); + } + + void init(const diameter_entry_t _simplex, const index_t _dim) { + idx_below = get_index(_simplex); + idx_above= 0; + v = parent.n - 1; + k = _dim + 1; + vertices.resize(_dim + 1); + simplex = _simplex; parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin()); } @@ -740,12 +861,12 @@ template <> class ripser<sparse_distance_matrix>::simplex_coboundary_enumerator const ripser& parent; 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; public: @@ -753,12 +874,23 @@ public: const ripser& _parent) : parent(_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) { + binomial_coeff(parent.binomial_coeff) { + 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()); + } + } + + void init(const diameter_entry_t _simplex, const index_t _dim) { + idx_below = get_index(_simplex); + idx_above= 0; + k = _dim + 1; + vertices.resize(_dim + 1); + simplex = _simplex; + parent.get_simplex_vertices(idx_below, _dim, parent.n, vertices.rbegin()); 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()); @@ -1096,14 +1228,13 @@ int main(int argc, char** argv) { max = -std::numeric_limits<value_t>::infinity(), max_finite = max; int num_edges = 0; + value_t enclosing_radius = std::numeric_limits<value_t>::infinity(); if (threshold == std::numeric_limits<value_t>::max()) { - value_t enclosing_radius = std::numeric_limits<value_t>::infinity(); for (size_t i = 0; i < dist.size(); ++i) { value_t r_i = -std::numeric_limits<value_t>::infinity(); for (size_t j = 0; j < dist.size(); ++j) r_i = std::max(r_i, dist(i, j)); enclosing_radius = std::min(enclosing_radius, r_i); } - threshold = enclosing_radius; } for (auto d : dist.distances) { @@ -1115,9 +1246,9 @@ int main(int argc, char** argv) { } std::cout << "value range: [" << min << "," << max_finite << "]" << std::endl; - if (threshold >= max) { - std::cout << "distance matrix with " << dist.size() << " points" << std::endl; - ripser<compressed_lower_distance_matrix>(std::move(dist), dim_max, threshold, ratio, + if (threshold == std::numeric_limits<value_t>::max()) { + std::cout << "distance matrix with " << dist.size() << " points, using threshold at enclosing radius " << enclosing_radius << std::endl; + ripser<compressed_lower_distance_matrix>(std::move(dist), dim_max, enclosing_radius, ratio, modulus) .compute_barcodes(); } else { |