diff options
| -rw-r--r-- | ripser.cpp | 163 |
1 files changed, 138 insertions, 25 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 @@ -279,6 +279,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(); } }; @@ -383,6 +391,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 { @@ -426,8 +435,92 @@ 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, entry_hash_map& pivot_column_index, index_t dim) { @@ -458,7 +551,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)}); } } @@ -500,8 +595,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()); @@ -546,7 +644,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()) { @@ -631,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)] % @@ -641,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 @@ -1090,14 +1202,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) { @@ -1109,10 +1220,12 @@ 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, - modulus) + 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 { std::cout << "sparse distance matrix with " << dist.size() << " points and " |