From d324d2a034b037e22b413bcf79f69014fd35bfc9 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 30 Jun 2019 00:29:04 +0200 Subject: entry hash map stores pivot coefficients --- ripser.cpp | 53 +++++++++++++++++++++++++++++++++-------------------- 1 file changed, 33 insertions(+), 20 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index fde30ad..3d90ea6 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -50,14 +50,16 @@ #ifdef USE_GOOGLE_HASHMAP #include -template class hash_map : public google::dense_hash_map { +template +class hash_map : public google::dense_hash_map { public: - explicit hash_map() : google::dense_hash_map() { this->set_empty_key(-1); } + explicit hash_map() : google::dense_hash_map() { this->set_empty_key(-1); } inline void reserve(size_t hint) { this->resize(hint); } }; #else -template class hash_map : public std::unordered_map {}; +template +class hash_map : public std::unordered_map {}; #endif typedef float value_t; @@ -129,6 +131,18 @@ std::ostream& operator<<(std::ostream& stream, const entry_t& e) { const entry_t& get_entry(const entry_t& e) { return e; } +struct entry_hash { + std::size_t operator()(const entry_t& e) const { return std::hash()(get_index(e)); } +}; + +struct equal_index { + bool operator()(const entry_t& e, const entry_t& f) const { + return get_index(e) == get_index(f); + } +}; + +typedef hash_map entry_hash_map; + typedef std::pair diameter_index_t; value_t get_diameter(const diameter_index_t& i) { return i.first; } index_t get_index(const diameter_index_t& i) { return i.second; } @@ -414,9 +428,9 @@ public: class simplex_coboundary_enumerator; - void assemble_columns_to_reduce(std::vector& simplices, - std::vector& columns_to_reduce, - hash_map& pivot_column_index, index_t dim) { + void assemble_columns_to_reduce( + std::vector& simplices, std::vector& columns_to_reduce, + entry_hash_map& pivot_column_index, index_t dim) { #ifdef INDICATE_PROGRESS std::cout << "\033[K" @@ -437,7 +451,7 @@ public: next_simplices.push_back(std::make_pair(get_diameter(coface), get_index(coface))); - if (pivot_column_index.find(get_index(coface)) == pivot_column_index.end()) + if (pivot_column_index.find(get_entry(coface)) == pivot_column_index.end()) columns_to_reduce.push_back( std::make_pair(get_diameter(coface), get_index(coface))); } @@ -492,9 +506,9 @@ public: } template - diameter_entry_t init_coboundary_and_get_pivot(diameter_entry_t simplex, - Column& working_coboundary, const index_t& dim, - hash_map& pivot_column_index) { + diameter_entry_t init_coboundary_and_get_pivot( + diameter_entry_t simplex, Column& working_coboundary, const index_t& dim, + entry_hash_map& pivot_column_index) { bool check_for_emergent_pair = true; coface_entries.clear(); simplex_coboundary_enumerator cofaces(simplex, dim, *this); @@ -503,7 +517,7 @@ public: if (get_diameter(coface) <= threshold) { coface_entries.push_back(coface); if (check_for_emergent_pair && (get_diameter(simplex) == get_diameter(coface))) { - if (pivot_column_index.find(get_index(coface)) == pivot_column_index.end()) + if (pivot_column_index.find(get_entry(coface)) == pivot_column_index.end()) return coface; check_for_emergent_pair = false; } @@ -532,7 +546,8 @@ public: } void compute_pairs(std::vector& columns_to_reduce, - hash_map& pivot_column_index, index_t dim) { + entry_hash_map& pivot_column_index, + index_t dim) { #ifdef PRINT_PERSISTENCE_PAIRS std::cout << "persistence intervals in dim " << dim << ":" << std::endl; @@ -565,10 +580,11 @@ public: while (true) { if (get_index(pivot) != -1) { - auto pair = pivot_column_index.find(get_index(pivot)); + auto pair = pivot_column_index.find(get_entry(pivot)); if (pair != pivot_column_index.end()) { index_t index_column_to_add = pair->second; - coefficient_t factor_column_to_add = modulus - get_coefficient(pivot); + entry_t other_pivot = pair->first; + coefficient_t factor_column_to_add = modulus - get_coefficient(pivot) * multiplicative_inverse[get_coefficient(other_pivot)] % modulus; add_coboundary(reduction_matrix.cbegin(index_column_to_add), reduction_matrix.cend(index_column_to_add), @@ -586,17 +602,14 @@ public: std::cout << " [" << diameter << "," << death << ")" << std::endl; } #endif - pivot_column_index.insert( - std::make_pair(get_index(pivot), index_column_to_reduce)); + pivot_column_index.insert(std::make_pair(get_entry(pivot), + index_column_to_reduce)); reduction_matrix.append_column(); - const coefficient_t inverse = - multiplicative_inverse[get_coefficient(pivot)]; while (true) { diameter_entry_t e = pop_pivot(working_reduction_column, modulus); if (get_index(e) == -1) break; - set_coefficient(e, inverse * get_coefficient(e) % modulus); assert(get_coefficient(e) > 0); reduction_matrix.push_back(e); } @@ -623,7 +636,7 @@ public: compute_dim_0_pairs(simplices, columns_to_reduce); for (index_t dim = 1; dim <= dim_max; ++dim) { - hash_map pivot_column_index; + entry_hash_map pivot_column_index; pivot_column_index.reserve(columns_to_reduce.size()); compute_pairs(columns_to_reduce, pivot_column_index, dim); -- cgit v1.2.3 From 4e02a1df4ed48bfabeab1d473d2afaf837855508 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 30 Jun 2019 01:08:19 +0200 Subject: don't store 1s on diagonal --- ripser.cpp | 66 ++++++++++++++++++++++++++++++++++++++------------------------ 1 file changed, 41 insertions(+), 25 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 3d90ea6..038f6d5 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -428,9 +428,9 @@ public: class simplex_coboundary_enumerator; - void assemble_columns_to_reduce( - std::vector& simplices, std::vector& columns_to_reduce, - entry_hash_map& pivot_column_index, index_t dim) { + void assemble_columns_to_reduce(std::vector& simplices, + std::vector& columns_to_reduce, + entry_hash_map& pivot_column_index, index_t dim) { #ifdef INDICATE_PROGRESS std::cout << "\033[K" @@ -506,9 +506,9 @@ public: } template - diameter_entry_t init_coboundary_and_get_pivot( - diameter_entry_t simplex, Column& working_coboundary, const index_t& dim, - entry_hash_map& pivot_column_index) { + diameter_entry_t init_coboundary_and_get_pivot(diameter_entry_t simplex, + Column& working_coboundary, const index_t& dim, + entry_hash_map& pivot_column_index) { bool check_for_emergent_pair = true; coface_entries.clear(); simplex_coboundary_enumerator cofaces(simplex, dim, *this); @@ -527,27 +527,39 @@ public: return get_pivot(working_coboundary, modulus); } - template - void add_coboundary(Iterator column_begin, Iterator column_end, - coefficient_t factor_column_to_add, Column& working_reduction_column, - Column& working_coboundary, const index_t& dim) { - coface_entries.clear(); + template + void add_coboundary(diameter_entry_t simplex, coefficient_t factor, Column& working_coboundary, + const index_t& dim) { + set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); + + simplex_coboundary_enumerator cofaces(simplex, dim, *this); + while (cofaces.has_next()) { + diameter_entry_t coface = cofaces.next(); + if (get_diameter(coface) <= threshold) working_coboundary.push(coface); + } + } + + template + void add_coboundary(compressed_sparse_matrix& reduction_matrix, + index_t index_column_to_add, coefficient_t factor, + Column& working_reduction_column, Column& working_coboundary, + const index_t& dim) { + auto column_begin = reduction_matrix.cbegin(index_column_to_add), + column_end = reduction_matrix.cend(index_column_to_add); for (auto it = column_begin; it != column_end; ++it) { diameter_entry_t simplex = *it; - set_coefficient(simplex, get_coefficient(simplex) * factor_column_to_add % modulus); + set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); working_reduction_column.push(simplex); simplex_coboundary_enumerator cofaces(simplex, dim, *this); while (cofaces.has_next()) { diameter_entry_t coface = cofaces.next(); - if (get_diameter(coface) <= threshold) coface_entries.push_back(coface); + if (get_diameter(coface) <= threshold) working_coboundary.push(coface); } } - for (auto coface : coface_entries) working_coboundary.push(coface); } void compute_pairs(std::vector& columns_to_reduce, - entry_hash_map& pivot_column_index, - index_t dim) { + entry_hash_map& pivot_column_index, index_t dim) { #ifdef PRINT_PERSISTENCE_PAIRS std::cout << "persistence intervals in dim " << dim << ":" << std::endl; @@ -573,7 +585,7 @@ public: greater_diameter_or_smaller_index> working_reduction_column, working_coboundary; - working_reduction_column.push(column_to_reduce); + // working_reduction_column.push(column_to_reduce); diameter_entry_t pivot = init_coboundary_and_get_pivot( column_to_reduce, working_coboundary, dim, pivot_column_index); @@ -582,14 +594,18 @@ public: if (get_index(pivot) != -1) { auto pair = pivot_column_index.find(get_entry(pivot)); if (pair != pivot_column_index.end()) { - index_t index_column_to_add = pair->second; entry_t other_pivot = pair->first; - coefficient_t factor_column_to_add = modulus - get_coefficient(pivot) * multiplicative_inverse[get_coefficient(other_pivot)] % modulus; + index_t index_column_to_add = pair->second; + diameter_entry_t column_to_add(columns_to_reduce[index_column_to_add], 1); + + coefficient_t factor = + modulus - get_coefficient(pivot) * + multiplicative_inverse[get_coefficient(other_pivot)] % + modulus; - add_coboundary(reduction_matrix.cbegin(index_column_to_add), - reduction_matrix.cend(index_column_to_add), - factor_column_to_add, working_reduction_column, - working_coboundary, dim); + add_coboundary(column_to_add, factor, working_coboundary, dim); + add_coboundary(reduction_matrix, index_column_to_add, factor, + working_reduction_column, working_coboundary, dim); pivot = get_pivot(working_coboundary, modulus); } else { @@ -602,8 +618,8 @@ public: std::cout << " [" << diameter << "," << death << ")" << std::endl; } #endif - pivot_column_index.insert(std::make_pair(get_entry(pivot), - index_column_to_reduce)); + pivot_column_index.insert( + std::make_pair(get_entry(pivot), index_column_to_reduce)); reduction_matrix.append_column(); -- cgit v1.2.3 From 9f5b862d8caf6b3ef48ae96b34d119b573452711 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 30 Jun 2019 23:11:43 +0200 Subject: don't store 1s on diagonal (fixed) --- ripser.cpp | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 038f6d5..a4b052a 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -528,10 +528,9 @@ public: } template - void add_coboundary(diameter_entry_t simplex, coefficient_t factor, Column& working_coboundary, - const index_t& dim) { - set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); - + void add_coboundary(diameter_entry_t simplex, Column& working_reduction_column, + Column& working_coboundary, const index_t& dim) { + working_reduction_column.push(simplex); simplex_coboundary_enumerator cofaces(simplex, dim, *this); while (cofaces.has_next()) { diameter_entry_t coface = cofaces.next(); @@ -544,17 +543,12 @@ public: index_t index_column_to_add, coefficient_t factor, Column& working_reduction_column, Column& working_coboundary, const index_t& dim) { - auto column_begin = reduction_matrix.cbegin(index_column_to_add), - column_end = reduction_matrix.cend(index_column_to_add); - for (auto it = column_begin; it != column_end; ++it) { + for (auto it = reduction_matrix.cbegin(index_column_to_add); + it != reduction_matrix.cend(index_column_to_add); ++it) { diameter_entry_t simplex = *it; set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); working_reduction_column.push(simplex); - simplex_coboundary_enumerator cofaces(simplex, dim, *this); - while (cofaces.has_next()) { - diameter_entry_t coface = cofaces.next(); - if (get_diameter(coface) <= threshold) working_coboundary.push(coface); - } + add_coboundary(simplex, working_reduction_column, working_coboundary, dim); } } @@ -596,14 +590,13 @@ public: if (pair != pivot_column_index.end()) { entry_t other_pivot = pair->first; index_t index_column_to_add = pair->second; - diameter_entry_t column_to_add(columns_to_reduce[index_column_to_add], 1); - coefficient_t factor = modulus - get_coefficient(pivot) * multiplicative_inverse[get_coefficient(other_pivot)] % modulus; + diameter_entry_t column_to_add(columns_to_reduce[index_column_to_add], factor); - add_coboundary(column_to_add, factor, working_coboundary, dim); + add_coboundary(column_to_add, working_reduction_column, working_coboundary, dim); add_coboundary(reduction_matrix, index_column_to_add, factor, working_reduction_column, working_coboundary, dim); -- cgit v1.2.3 From d9ff05f213a63add6d7cb43f69e63b6d84905177 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Tue, 2 Jul 2019 00:18:20 +0200 Subject: added some const specifiers --- ripser.cpp | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index a4b052a..65d84d7 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -247,7 +247,7 @@ struct sparse_distance_matrix { : neighbors(std::move(_neighbors)), num_edges(_num_edges) {} template - sparse_distance_matrix(const DistanceMatrix& mat, value_t threshold) + sparse_distance_matrix(const DistanceMatrix& mat, const value_t threshold) : neighbors(mat.size()), num_edges(0) { for (index_t i = 0; i < size(); ++i) @@ -285,7 +285,7 @@ class union_find { std::vector rank; public: - union_find(index_t n) : parent(n), rank(n, 0) { + union_find(const index_t n) : parent(n), rank(n, 0) { for (index_t i = 0; i < n; ++i) parent[i] = i; } @@ -326,7 +326,7 @@ template diameter_entry_t pop_pivot(Column& column, coefficien return pivot; } -template diameter_entry_t get_pivot(Column& column, coefficient_t modulus) { +template diameter_entry_t get_pivot(Column& column, const coefficient_t modulus) { diameter_entry_t result = pop_pivot(column, modulus); if (get_index(result) != -1) column.push(result); return result; @@ -344,7 +344,7 @@ public: return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1]; } - typename std::vector::const_iterator cend(size_t index) const { + typename std::vector::const_iterator cend(const size_t index) const { assert(index < size()); return entries.cbegin() + bounds[index]; } @@ -356,7 +356,7 @@ public: void append_column() { bounds.push_back(entries.size()); } - void push_back(ValueType e) { + void push_back(const ValueType e) { assert(0 < size()); entries.push_back(e); ++bounds.back(); @@ -373,7 +373,8 @@ public: } }; -template index_t get_max(index_t top, index_t bottom, const Predicate pred) { +template +index_t get_max(index_t top, const index_t bottom, const Predicate pred) { if (!pred(top)) { index_t count = top - bottom; while (count > 0) { @@ -506,7 +507,7 @@ public: } template - diameter_entry_t init_coboundary_and_get_pivot(diameter_entry_t simplex, + 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; @@ -528,8 +529,8 @@ public: } template - void add_coboundary(diameter_entry_t simplex, Column& working_reduction_column, - Column& working_coboundary, const index_t& dim) { + void add_coboundary(const diameter_entry_t simplex, Column& working_reduction_column, + Column& working_coboundary, const index_t& dim) { working_reduction_column.push(simplex); simplex_coboundary_enumerator cofaces(simplex, dim, *this); while (cofaces.has_next()) { @@ -540,7 +541,7 @@ public: template void add_coboundary(compressed_sparse_matrix& reduction_matrix, - index_t index_column_to_add, coefficient_t factor, + const index_t index_column_to_add, const coefficient_t factor, Column& working_reduction_column, Column& working_coboundary, const index_t& dim) { for (auto it = reduction_matrix.cbegin(index_column_to_add); @@ -553,7 +554,7 @@ public: } void compute_pairs(std::vector& columns_to_reduce, - entry_hash_map& pivot_column_index, index_t dim) { + entry_hash_map& pivot_column_index, const index_t dim) { #ifdef PRINT_PERSISTENCE_PAIRS std::cout << "persistence intervals in dim " << dim << ":" << std::endl; @@ -594,9 +595,11 @@ public: modulus - get_coefficient(pivot) * multiplicative_inverse[get_coefficient(other_pivot)] % modulus; - diameter_entry_t column_to_add(columns_to_reduce[index_column_to_add], factor); + diameter_entry_t column_to_add(columns_to_reduce[index_column_to_add], + factor); - add_coboundary(column_to_add, working_reduction_column, working_coboundary, dim); + add_coboundary(column_to_add, working_reduction_column, working_coboundary, + dim); add_coboundary(reduction_matrix, index_column_to_add, factor, working_reduction_column, working_coboundary, dim); @@ -709,7 +712,7 @@ template <> class ripser::simplex_coboundary_enumerator index_diameter_t neighbor; public: - simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, + simplex_coboundary_enumerator(const diameter_entry_t _simplex, const index_t _dim, const ripser& _parent) : parent(_parent), idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), k(_dim + 1), max_vertex_below(parent.n - 1), simplex(_simplex), modulus(parent.modulus), @@ -919,7 +922,7 @@ compressed_lower_distance_matrix read_binary(std::istream& input_stream) { return compressed_lower_distance_matrix(std::move(distances)); } -compressed_lower_distance_matrix read_file(std::istream& input_stream, file_format format) { +compressed_lower_distance_matrix read_file(std::istream& input_stream, const file_format format) { switch (format) { case LOWER_DISTANCE_MATRIX: return read_lower_distance_matrix(input_stream); -- cgit v1.2.3 From 4e2129b11a4dc491f218bc9840d1f19eb898a391 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Tue, 2 Jul 2019 11:06:16 +0200 Subject: fixed google hash map support --- ripser.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/ripser.cpp b/ripser.cpp index 65d84d7..55e94b5 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -108,6 +108,7 @@ struct __attribute__((packed)) entry_t { coefficient_t coefficient; entry_t(index_t _index, coefficient_t _coefficient) : index(_index), coefficient(_coefficient) {} + entry_t(index_t _index) : index(_index), coefficient(0) {} entry_t() : index(0), coefficient(0) {} }; -- cgit v1.2.3