From 246f1fe3efcb794b88c95865477be2b7dc0cf709 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 20 Jun 2019 01:26:50 +0200 Subject: consolidation --- ripser.cpp | 253 +++++++++++++++++++++++-------------------------------------- 1 file changed, 94 insertions(+), 159 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 37b4987..3d08300 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -205,6 +205,9 @@ public: index_t num_edges; + mutable std::vector::const_reverse_iterator> neighbor_it; + mutable std::vector::const_reverse_iterator> neighbor_end; + sparse_distance_matrix(std::vector>&& _neighbors, index_t _num_edges) : neighbors(std::move(_neighbors)), num_edges(_num_edges) {} @@ -391,8 +394,6 @@ template class ripser { const binomial_coeff_table binomial_coeff; std::vector multiplicative_inverse; mutable std::vector vertices; - mutable std::vector::const_reverse_iterator> neighbor_it; - mutable std::vector::const_reverse_iterator> neighbor_end; mutable std::vector coface_entries; public: @@ -424,24 +425,49 @@ public: return out; } - value_t compute_diameter(const index_t index, index_t dim) const { - value_t diam = -std::numeric_limits::infinity(); + 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) { + +#ifdef INDICATE_PROGRESS + std::cout << "\033[K" + << "assembling columns" << std::flush << "\r"; +#endif + + --dim; + columns_to_reduce.clear(); + + std::vector next_simplices; + + for (diameter_index_t simplex : simplices) { + simplex_coboundary_enumerator cofaces(diameter_entry_t(simplex, 1), dim, *this); - vertices.clear(); - get_simplex_vertices(index, dim, dist.size(), std::back_inserter(vertices)); + while (cofaces.has_next(false)) { + auto coface = cofaces.next(); - 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])); + 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()) + columns_to_reduce.push_back( + std::make_pair(get_diameter(coface), get_index(coface))); } - return diam; - } + } - class simplex_coboundary_enumerator; + simplices.swap(next_simplices); - void assemble_columns_to_reduce(std::vector& simplices, - std::vector& columns_to_reduce, - hash_map& pivot_column_index, index_t dim); +#ifdef INDICATE_PROGRESS + std::cout << "\033[K" + << "sorting " << columns_to_reduce.size() << " columns" << std::flush << "\r"; +#endif + + std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), + greater_diameter_or_smaller_index()); +#ifdef INDICATE_PROGRESS + std::cout << "\033[K"; +#endif + } void compute_dim_0_pairs(std::vector& edges, std::vector& columns_to_reduce) { @@ -458,8 +484,7 @@ public: std::vector vertices_of_edge(2); for (auto e : edges) { - vertices_of_edge.clear(); - get_simplex_vertices(get_index(e), 1, n, std::back_inserter(vertices_of_edge)); + get_simplex_vertices(get_index(e), 1, n, vertices_of_edge.rbegin()); index_t u = dset.find(vertices_of_edge[0]), v = dset.find(vertices_of_edge[1]); if (u != v) { @@ -485,7 +510,33 @@ public: Column& working_reduction_column, Column& working_coboundary, const index_t& dim, hash_map& pivot_column_index, - bool& might_be_apparent_pair); + bool& might_be_apparent_pair) { + 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); + + working_reduction_column.push(simplex); + + coface_entries.clear(); + 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 (might_be_apparent_pair && (get_diameter(simplex) == get_diameter(coface))) { + if (pivot_column_index.find(get_index(coface)) == + pivot_column_index.end()) { + return coface; + } + might_be_apparent_pair = false; + } + } + } + for (auto coface : coface_entries) working_coboundary.push(coface); + } + + return get_pivot(working_coboundary, modulus); + } void compute_pairs(std::vector& columns_to_reduce, hash_map& pivot_column_index, index_t dim) { @@ -533,11 +584,11 @@ public: bool might_be_apparent_pair = (index_column_to_reduce == index_column_to_add); while (true) { - pivot = add_coboundary_and_get_pivot( - reduction_matrix.cbegin(index_column_to_add), - reduction_matrix.cend(index_column_to_add), - factor_column_to_add, working_reduction_column, - working_coboundary, dim, pivot_column_index, might_be_apparent_pair); + pivot = add_coboundary_and_get_pivot(reduction_matrix.cbegin(index_column_to_add), + reduction_matrix.cend(index_column_to_add), + factor_column_to_add, working_reduction_column, + working_coboundary, dim, pivot_column_index, + might_be_apparent_pair); if (get_index(pivot) != -1) { auto pair = pivot_column_index.find(get_index(pivot)); @@ -627,11 +678,12 @@ public: : 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.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.begin()); + parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin()); } - bool has_next() { + bool has_next(bool all_cofaces = true) { while ((v != -1) && (binomial_coeff(v, k) <= idx_below)) { + if (!all_cofaces) return false; idx_below -= binomial_coeff(v, k); idx_above += binomial_coeff(v, k + 1); --v; @@ -651,38 +703,6 @@ public: } }; -template -template -diameter_entry_t ripser::add_coboundary_and_get_pivot( - Iterator column_begin, Iterator column_end, coefficient_t factor_column_to_add, - Column& working_reduction_column, Column& working_coboundary, const index_t& dim, - hash_map& pivot_column_index, bool& might_be_apparent_pair) { - 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); - - working_reduction_column.push(simplex); - - coface_entries.clear(); - 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 (might_be_apparent_pair && (get_diameter(simplex) == get_diameter(coface))) { - if (pivot_column_index.find(get_index(coface)) == pivot_column_index.end()) { - return coface; - } - might_be_apparent_pair = false; - } - } - } - for (auto coface : coface_entries) working_coboundary.push(coface); - } - - return get_pivot(working_coboundary, modulus); -} - template <> class ripser::simplex_coboundary_enumerator { private: const ripser& parent; @@ -693,7 +713,7 @@ private: const sparse_distance_matrix& dist; const binomial_coeff_table& binomial_coeff; - std::vector& vertices; + std::vector vertices; std::vector::const_reverse_iterator>& neighbor_it; std::vector::const_reverse_iterator>& neighbor_end; index_diameter_t x; @@ -703,14 +723,13 @@ public: 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), - dist(parent.dist), binomial_coeff(parent.binomial_coeff), vertices(parent.vertices), - neighbor_it(parent.neighbor_it), neighbor_end(parent.neighbor_end) { + dist(parent.dist), binomial_coeff(parent.binomial_coeff), vertices(_dim + 1), + neighbor_it(dist.neighbor_it), neighbor_end(dist.neighbor_end) { neighbor_it.clear(); neighbor_end.clear(); - vertices.clear(); - parent.get_simplex_vertices(idx_below, _dim, parent.n, std::back_inserter(vertices)); + parent.get_simplex_vertices(idx_below, _dim, parent.n, vertices.rbegin()); for (auto v : vertices) { neighbor_it.push_back(dist.neighbors[v].rbegin()); @@ -731,8 +750,13 @@ public: else x = std::max(x, y); } - return all_cofaces || !(k > 0 && parent.get_next_vertex(max_vertex_below, idx_below, - k) > get_index(x)); + while (k > 0 && vertices[k - 1] > get_index(x)) { + if (!all_cofaces) return false; + idx_below -= binomial_coeff(vertices[k - 1], k); + idx_above += binomial_coeff(vertices[k - 1], k + 1); + --k; + } + return true; continue_outer:; } return false; @@ -740,29 +764,21 @@ public: diameter_entry_t next() { ++neighbor_it[0]; - - while (k > 0 && parent.get_next_vertex(max_vertex_below, idx_below, k) > get_index(x)) { - idx_below -= binomial_coeff(max_vertex_below, k); - idx_above += binomial_coeff(max_vertex_below, k + 1); - --k; - } - value_t coface_diameter = std::max(get_diameter(simplex), get_diameter(x)); - + index_t coface_index = idx_above + binomial_coeff(get_index(x), k + 1) + idx_below; coefficient_t coface_coefficient = (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; - - return diameter_entry_t(coface_diameter, - idx_above + binomial_coeff(get_index(x), k + 1) + idx_below, - coface_coefficient); + return diameter_entry_t(coface_diameter, coface_index, coface_coefficient); } }; template <> std::vector ripser::get_edges() { std::vector edges; + std::vector vertices(2); for (index_t index = binomial_coeff(n, 2); index-- > 0;) { - value_t diameter = compute_diameter(index, 1); - if (diameter <= threshold) edges.push_back(std::make_pair(diameter, index)); + get_simplex_vertices(index, 1, dist.size(), vertices.rbegin()); + value_t length = dist(vertices[0], vertices[1]); + if (length <= threshold) edges.push_back(std::make_pair(length, index)); } return edges; } @@ -777,87 +793,6 @@ template <> std::vector ripser::get_ed return edges; } -template <> -void ripser::assemble_columns_to_reduce( - std::vector& simplices, std::vector& columns_to_reduce, - hash_map& pivot_column_index, index_t dim) { - index_t num_simplices = binomial_coeff(n, dim + 1); - - columns_to_reduce.clear(); - -#ifdef INDICATE_PROGRESS - std::cout << "\033[K" - << "assembling " << num_simplices << " columns" << std::flush << "\r"; -#endif - - for (index_t index = 0; index < num_simplices; ++index) { - if (pivot_column_index.find(index) == pivot_column_index.end()) { - value_t diameter = compute_diameter(index, dim); - if (diameter <= threshold) columns_to_reduce.push_back(std::make_pair(diameter, index)); -#ifdef INDICATE_PROGRESS - if ((index + 1) % 1000000 == 0) - std::cout << "\033[K" - << "assembled " << columns_to_reduce.size() << " out of " << (index + 1) - << "/" << num_simplices << " columns" << std::flush << "\r"; -#endif - } - } - -#ifdef INDICATE_PROGRESS - std::cout << "\033[K" - << "sorting " << num_simplices << " columns" << std::flush << "\r"; -#endif - - std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), - greater_diameter_or_smaller_index()); -#ifdef INDICATE_PROGRESS - std::cout << "\033[K"; -#endif -} - -template <> -void ripser::assemble_columns_to_reduce( - std::vector& simplices, std::vector& columns_to_reduce, - hash_map& pivot_column_index, index_t dim) { - -#ifdef INDICATE_PROGRESS - std::cout << "\033[K" - << "assembling columns" << std::flush << "\r"; -#endif - - --dim; - columns_to_reduce.clear(); - - std::vector next_simplices; - - for (diameter_index_t simplex : simplices) { - simplex_coboundary_enumerator cofaces(diameter_entry_t(simplex, 1), dim, *this); - - while (cofaces.has_next(false)) { - auto coface = cofaces.next(); - - 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()) - columns_to_reduce.push_back( - std::make_pair(get_diameter(coface), get_index(coface))); - } - } - - simplices.swap(next_simplices); - -#ifdef INDICATE_PROGRESS - std::cout << "\033[K" - << "sorting " << columns_to_reduce.size() << " columns" << std::flush << "\r"; -#endif - - std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), - greater_diameter_or_smaller_index()); -#ifdef INDICATE_PROGRESS - std::cout << "\033[K"; -#endif -} - enum file_format { LOWER_DISTANCE_MATRIX, UPPER_DISTANCE_MATRIX, -- cgit v1.2.3 From cb6163df1dd96e2f65c5b4506fd88c1266562e86 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Mon, 24 Jun 2019 19:24:34 +0200 Subject: structs instead of classes --- README.md | 21 +++++++++++---------- ripser.cpp | 62 ++++++++++++++++++++++---------------------------------------- 2 files changed, 33 insertions(+), 50 deletions(-) diff --git a/README.md b/README.md index c8f82cc..3c556c0 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # Ripser -Copyright © 2015–2018 [Ulrich Bauer]. +Copyright © 2015–2019 [Ulrich Bauer]. ### Description @@ -24,6 +24,8 @@ Input formats currently supported by Ripser: - comma-separated values upper triangular distance matrix (MATLAB output from the function `pdist`) - comma-separated values full distance matrix - [DIPHA] distance matrix data + - sparse distance matrix in Sparse Triplet format + - binary lower triangular distance matrix - point cloud data Ripser's efficiency is based on a few important concepts and principles, building on key previous and concurrent developments by other researchers in computational topology: @@ -37,7 +39,7 @@ Ripser's efficiency is based on a few important concepts and principles, buildin ### Version -[Latest release][latest-release]: 1.0.1 (September 2016) +[Latest release][latest-release]: 1.1 (June 2019) ### Building @@ -56,16 +58,14 @@ make Ripser supports several compile-time options. They are switched on by defining the C preprocessor macros listed below, either using `#define` in the code or by passing an argument to the compiler. The following options are supported: - - `ASSEMBLE_REDUCTION_MATRIX`: store the reduction matrix; may affect computation time but also memory usage; recommended for large and difficult problem instances - - `USE_COEFFICIENTS`: enable support for coefficients in a prime field - `INDICATE_PROGRESS`: indicate the current progress in the console - `PRINT_PERSISTENCE_PAIRS`: output the computed persistence pairs (enabled by default in the code; comment out to disable) - `USE_GOOGLE_HASHMAP`: enable support for Google's [sparsehash] data structure; may further reduce memory footprint -For example, to build Ripser with support for coefficients: +For example, to build Ripser with support for Google's hashmap: ```sh -$ c++ -std=c++11 ripser.cpp -o ripser -Ofast -D NDEBUG -D USE_COEFFICIENTS +$ c++ -std=c++11 ripser.cpp -o ripser -Ofast -D NDEBUG -D USE_GOOGLE_HASHMAP ``` A Makefile is provided with some variants of the above options. Use `make all` to build them. The default `make` builds a binary with the default options. @@ -78,19 +78,20 @@ The input is given either in a file whose name is passed as an argument, or thro - `distance`: full distance matrix; similar to the above, but for all entries of the distance matrix. One line per row of the matrix; only the part below the diagonal is actually read. - `dipha`: DIPHA distance matrix as described on the [DIPHA] website. - `point-cloud`: point cloud; a comma (or whitespace, or other non-numerical character) separated list of coordinates of the points in some Euclidean space, one point per line. + - `sparse`: sparse distance matrix in Sparse Triplet format + - `binary`: lower distance matrix in binary file format; a sequence of the distance matrix entries below the diagonal in 64 bit double format (IEEE 754, little endian). - `--dim k`: compute persistent homology up to dimension *k*. - `--threshold t`: compute Rips complexes up to diameter *t*. - - `--modulus p`: compute homology with coefficients in the prime field Z/*p*Z (only available when built with the option `USE_COEFFICIENTS`). + - `--modulus p`: compute homology with coefficients in the prime field Z/*p*Z. + - `--ratio `: only show persistence pairs with death/birth ratio > *r*. - -### Planned features +### Experimental features The following features are currently planned for future versions: - computation of representative cycles for persistent homology (currenly only *co*cycles are computed) - - support for sparse distance matrices Prototype implementations are already avaliable; please contact the author if one of these features might be relevant for your research. diff --git a/ripser.cpp b/ripser.cpp index 3d08300..bdb452a 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -41,9 +41,6 @@ //#define USE_GOOGLE_HASHMAP -#include -#include -#include #include #include #include @@ -141,8 +138,7 @@ typedef std::pair index_diameter_t; index_t get_index(const index_diameter_t& i) { return i.first; } value_t get_diameter(const index_diameter_t& i) { return i.second; } -class diameter_entry_t : public std::pair { -public: +struct diameter_entry_t : std::pair { diameter_entry_t() {} diameter_entry_t(value_t _diameter, index_t _index, coefficient_t _coefficient) : std::pair(_diameter, make_entry(_index, _coefficient)) {} @@ -172,8 +168,7 @@ template struct greater_diameter_or_smaller_index { enum compressed_matrix_layout { LOWER_TRIANGULAR, UPPER_TRIANGULAR }; -template class compressed_distance_matrix { -public: +template struct compressed_distance_matrix { std::vector distances; std::vector rows; @@ -199,8 +194,7 @@ public: size_t size() const { return rows.size(); } }; -class sparse_distance_matrix { -public: +struct sparse_distance_matrix { std::vector> neighbors; index_t num_edges; @@ -258,8 +252,7 @@ value_t compressed_distance_matrix::operator()(const index_t i typedef compressed_distance_matrix compressed_lower_distance_matrix; typedef compressed_distance_matrix compressed_upper_distance_matrix; -class euclidean_distance_matrix { -public: +struct euclidean_distance_matrix { std::vector> points; euclidean_distance_matrix(std::vector>&& _points) @@ -393,7 +386,6 @@ template class ripser { coefficient_t modulus; const binomial_coeff_table binomial_coeff; std::vector multiplicative_inverse; - mutable std::vector vertices; mutable std::vector coface_entries; public: @@ -474,7 +466,6 @@ public: union_find dset(n); edges = get_edges(); - std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index()); @@ -553,11 +544,7 @@ public: std::priority_queue, greater_diameter_or_smaller_index> - working_reduction_column; - - std::priority_queue, - greater_diameter_or_smaller_index> - working_coboundary; + working_reduction_column, working_coboundary; value_t diameter = get_diameter(column_to_reduce); @@ -644,7 +631,6 @@ public: std::vector get_edges(); void compute_barcodes() { - std::vector simplices, columns_to_reduce; compute_dim_0_pairs(simplices, columns_to_reduce); @@ -655,16 +641,14 @@ public: compute_pairs(columns_to_reduce, pivot_column_index, dim); - if (dim < dim_max) { + if (dim < dim_max) assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, dim + 1); - } } } }; template <> class ripser::simplex_coboundary_enumerator { -private: index_t idx_below, idx_above, v, k; std::vector vertices; const diameter_entry_t simplex; @@ -704,16 +688,15 @@ public: }; template <> class ripser::simplex_coboundary_enumerator { -private: const ripser& parent; index_t idx_below, idx_above, v, k, max_vertex_below; + std::vector vertices; const diameter_entry_t simplex; const coefficient_t modulus; const sparse_distance_matrix& dist; const binomial_coeff_table& binomial_coeff; - std::vector vertices; std::vector::const_reverse_iterator>& neighbor_it; std::vector::const_reverse_iterator>& neighbor_end; index_diameter_t x; @@ -800,7 +783,7 @@ enum file_format { POINT_CLOUD, DIPHA, SPARSE, - RIPSER + BINARY }; template T read(std::istream& s) { @@ -841,9 +824,7 @@ compressed_lower_distance_matrix read_point_cloud(std::istream& input_stream) { } sparse_distance_matrix read_sparse_distance_matrix(std::istream& input_stream) { - std::vector> neighbors; - index_t num_edges = 0; std::string line; @@ -931,7 +912,7 @@ compressed_lower_distance_matrix read_dipha(std::istream& input_stream) { return compressed_lower_distance_matrix(std::move(distances)); } -compressed_lower_distance_matrix read_ripser(std::istream& input_stream) { +compressed_lower_distance_matrix read_binary(std::istream& input_stream) { std::vector distances; while (!input_stream.eof()) distances.push_back(read(input_stream)); return compressed_lower_distance_matrix(std::move(distances)); @@ -950,7 +931,7 @@ compressed_lower_distance_matrix read_file(std::istream& input_stream, file_form case DIPHA: return read_dipha(input_stream); default: - return read_ripser(input_stream); + return read_binary(input_stream); } } @@ -973,18 +954,19 @@ void print_usage_and_exit(int exit_code) { << " dipha (distance matrix in DIPHA file format)" << std::endl << " sparse (sparse distance matrix in Sparse Triplet format)" << std::endl - << " ripser (distance matrix in Ripser binary file format)" + << " binary (lower triangular distance matrix in binary format)" << std::endl - << " --dim compute persistent homology up to dimension " << std::endl - << " --threshold compute Rips complexes up to diameter " << std::endl - << " --modulus

compute homology with coefficients in the prime field Z/

Z" + << " --dim compute persistent homology up to dimension k" << std::endl + << " --threshold compute Rips complexes up to diameter t" << std::endl + << " --modulus

compute homology with coefficients in the prime field Z/pZ" + << std::endl + << " --ratio only show persistence pairs with death/birth ratio > r" << std::endl; exit(exit_code); } int main(int argc, char** argv) { - const char* filename = nullptr; file_format format = DISTANCE_MATRIX; @@ -1017,20 +999,20 @@ int main(int argc, char** argv) { if (next_pos != parameter.size()) print_usage_and_exit(-1); } else if (arg == "--format") { std::string parameter = std::string(argv[++i]); - if (parameter == "lower-distance") + if (parameter.rfind("lower", 0) == 0) format = LOWER_DISTANCE_MATRIX; - else if (parameter == "upper-distance") + else if (parameter.rfind("upper", 0) == 0) format = UPPER_DISTANCE_MATRIX; - else if (parameter == "distance") + else if (parameter.rfind("dist", 0) == 0) format = DISTANCE_MATRIX; - else if (parameter == "point-cloud") + else if (parameter.rfind("point", 0) == 0) format = POINT_CLOUD; else if (parameter == "dipha") format = DIPHA; else if (parameter == "sparse") format = SPARSE; - else if (parameter == "ripser") - format = RIPSER; + else if (parameter == "binary") + format = BINARY; else print_usage_and_exit(-1); } else if (arg == "--modulus") { -- cgit v1.2.3 From 4ece0bfaaab2affee63ec8dd5e55345ddb73e3c2 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Tue, 25 Jun 2019 08:54:30 +0200 Subject: cleanup --- ripser.cpp | 91 ++++++++++++++++++++++---------------------------------------- 1 file changed, 32 insertions(+), 59 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index bdb452a..2455291 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -172,8 +172,6 @@ template struct compressed_distance_matrix { std::vector distances; std::vector rows; - void init_rows(); - compressed_distance_matrix(std::vector&& _distances) : distances(std::move(_distances)), rows((1 + std::sqrt(1 + 8 * distances.size())) / 2) { assert(distances.size() == size() * (size() - 1) / 2); @@ -190,8 +188,8 @@ template struct compressed_distance_matrix { } value_t operator()(const index_t i, const index_t j) const; - size_t size() const { return rows.size(); } + void init_rows(); }; struct sparse_distance_matrix { @@ -289,6 +287,7 @@ public: } return z; } + void link(index_t x, index_t y) { if ((x = find(x)) == (y = find(y))) return; if (rank[x] > rank[y]) @@ -425,7 +424,8 @@ public: #ifdef INDICATE_PROGRESS std::cout << "\033[K" - << "assembling columns" << std::flush << "\r"; + << "assembling columns" + << "\r" << std::flush; #endif --dim; @@ -451,7 +451,8 @@ public: #ifdef INDICATE_PROGRESS std::cout << "\033[K" - << "sorting " << columns_to_reduce.size() << " columns" << std::flush << "\r"; + << "sorting " << columns_to_reduce.size() << " columns" + << "\r" << std::flush; #endif std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), @@ -463,16 +464,15 @@ public: void compute_dim_0_pairs(std::vector& edges, std::vector& columns_to_reduce) { +#ifdef PRINT_PERSISTENCE_PAIRS + std::cout << "persistence intervals in dim 0:" << std::endl; +#endif + union_find dset(n); edges = get_edges(); std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index()); - -#ifdef PRINT_PERSISTENCE_PAIRS - std::cout << "persistence intervals in dim 0:" << std::endl; -#endif - std::vector vertices_of_edge(2); for (auto e : edges) { get_simplex_vertices(get_index(e), 1, n, vertices_of_edge.rbegin()); @@ -491,7 +491,7 @@ public: #ifdef PRINT_PERSISTENCE_PAIRS for (index_t i = 0; i < n; ++i) - if (dset.find(i) == i) std::cout << " [0, )" << std::endl << std::flush; + if (dset.find(i) == i) std::cout << " [0, )" << std::endl; #endif } @@ -515,17 +515,14 @@ public: if (get_diameter(coface) <= threshold) { coface_entries.push_back(coface); if (might_be_apparent_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_index(coface)) == pivot_column_index.end()) return coface; - } might_be_apparent_pair = false; } } } for (auto coface : coface_entries) working_coboundary.push(coface); } - return get_pivot(working_coboundary, modulus); } @@ -541,11 +538,6 @@ public: for (index_t index_column_to_reduce = 0; index_column_to_reduce < columns_to_reduce.size(); ++index_column_to_reduce) { auto column_to_reduce = columns_to_reduce[index_column_to_reduce]; - - std::priority_queue, - greater_diameter_or_smaller_index> - working_reduction_column, working_coboundary; - value_t diameter = get_diameter(column_to_reduce); #ifdef INDICATE_PROGRESS @@ -553,13 +545,11 @@ public: std::cout << "\033[K" << "reducing column " << index_column_to_reduce + 1 << "/" << columns_to_reduce.size() << " (diameter " << diameter << ")" - << std::flush << "\r"; + << "\r" << std::flush; #endif index_t index_column_to_add = index_column_to_reduce; - diameter_entry_t pivot; - // start with factor 1 in order to initialize working_coboundary // with the coboundary of the simplex with index column_to_reduce coefficient_t factor_column_to_add = 1; @@ -568,18 +558,19 @@ public: reduction_matrix.append_column(); reduction_matrix.push_back(diameter_entry_t(column_to_reduce, 1)); + std::priority_queue, + greater_diameter_or_smaller_index> + working_reduction_column, working_coboundary; bool might_be_apparent_pair = (index_column_to_reduce == index_column_to_add); - + diameter_entry_t pivot; while (true) { pivot = add_coboundary_and_get_pivot(reduction_matrix.cbegin(index_column_to_add), reduction_matrix.cend(index_column_to_add), factor_column_to_add, working_reduction_column, working_coboundary, dim, pivot_column_index, might_be_apparent_pair); - if (get_index(pivot) != -1) { auto pair = pivot_column_index.find(get_index(pivot)); - if (pair != pivot_column_index.end()) { index_column_to_add = pair->second; factor_column_to_add = modulus - get_coefficient(pivot); @@ -590,21 +581,19 @@ public: #ifdef INDICATE_PROGRESS std::cout << "\033[K"; #endif - std::cout << " [" << diameter << "," << death << ")" << std::endl - << std::flush; + std::cout << " [" << diameter << "," << death << ")" << std::endl; } #endif pivot_column_index.insert( std::make_pair(get_index(pivot), index_column_to_reduce)); - const coefficient_t inverse = - multiplicative_inverse[get_coefficient(pivot)]; - // replace current column of reduction_matrix (with a single diagonal 1 // entry) by reduction_column (possibly with a different entry on the // diagonal) reduction_matrix.pop_back(); + 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; @@ -616,13 +605,12 @@ public: } } else { #ifdef PRINT_PERSISTENCE_PAIRS - std::cout << " [" << diameter << ", )" << std::endl << std::flush; + std::cout << " [" << diameter << ", )" << std::endl; #endif break; } } } - #ifdef INDICATE_PROGRESS std::cout << "\033[K"; #endif @@ -689,17 +677,15 @@ public: template <> class ripser::simplex_coboundary_enumerator { const ripser& parent; - index_t idx_below, idx_above, v, k, max_vertex_below; std::vector vertices; const diameter_entry_t simplex; const coefficient_t modulus; const sparse_distance_matrix& dist; const binomial_coeff_table& binomial_coeff; - std::vector::const_reverse_iterator>& neighbor_it; std::vector::const_reverse_iterator>& neighbor_end; - index_diameter_t x; + index_diameter_t neighbor; public: simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, @@ -708,12 +694,10 @@ public: k(_dim + 1), max_vertex_below(parent.n - 1), simplex(_simplex), modulus(parent.modulus), dist(parent.dist), binomial_coeff(parent.binomial_coeff), vertices(_dim + 1), 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()); @@ -722,18 +706,17 @@ public: bool has_next(bool all_cofaces = true) { for (auto &it0 = neighbor_it[0], &end0 = neighbor_end[0]; it0 != end0; ++it0) { - x = *it0; + neighbor = *it0; for (size_t idx = 1; idx < neighbor_it.size(); ++idx) { auto &it = neighbor_it[idx], end = neighbor_end[idx]; - while (get_index(*it) > get_index(x)) + while (get_index(*it) > get_index(neighbor)) if (++it == end) return false; - auto y = *it; - if (get_index(y) != get_index(x)) + if (get_index(*it) != get_index(neighbor)) goto continue_outer; else - x = std::max(x, y); + neighbor = std::max(neighbor, *it); } - while (k > 0 && vertices[k - 1] > get_index(x)) { + while (k > 0 && vertices[k - 1] > get_index(neighbor)) { if (!all_cofaces) return false; idx_below -= binomial_coeff(vertices[k - 1], k); idx_above += binomial_coeff(vertices[k - 1], k + 1); @@ -747,8 +730,8 @@ public: diameter_entry_t next() { ++neighbor_it[0]; - value_t coface_diameter = std::max(get_diameter(simplex), get_diameter(x)); - index_t coface_index = idx_above + binomial_coeff(get_index(x), k + 1) + idx_below; + value_t coface_diameter = std::max(get_diameter(simplex), get_diameter(neighbor)); + index_t coface_index = idx_above + binomial_coeff(get_index(neighbor), k + 1) + idx_below; coefficient_t coface_coefficient = (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; return diameter_entry_t(coface_diameter, coface_index, coface_coefficient); @@ -786,9 +769,9 @@ enum file_format { BINARY }; -template T read(std::istream& s) { +template T read(std::istream& input_stream) { T result; - s.read(reinterpret_cast(&result), sizeof(T)); + input_stream.read(reinterpret_cast(&result), sizeof(T)); return result; // on little endian: boost::endian::little_to_native(result); } @@ -809,14 +792,11 @@ compressed_lower_distance_matrix read_point_cloud(std::istream& input_stream) { } euclidean_distance_matrix eucl_dist(std::move(points)); - index_t n = eucl_dist.size(); - std::cout << "point cloud with " << n << " points in dimension " << eucl_dist.points.front().size() << std::endl; std::vector distances; - for (int i = 0; i < n; ++i) for (int j = 0; j < i; ++j) distances.push_back(eucl_dist(i, j)); @@ -960,9 +940,7 @@ void print_usage_and_exit(int exit_code) { << " --threshold compute Rips complexes up to diameter t" << std::endl << " --modulus

compute homology with coefficients in the prime field Z/pZ" << std::endl - << " --ratio only show persistence pairs with death/birth ratio > r" - << std::endl; - + << " --ratio only show persistence pairs with death/birth ratio > r" << std::endl; exit(exit_code); } @@ -973,9 +951,7 @@ int main(int argc, char** argv) { index_t dim_max = 1; value_t threshold = std::numeric_limits::max(); - float ratio = 1; - coefficient_t modulus = 2; for (index_t i = 1; i < argc; ++i) { @@ -1042,7 +1018,6 @@ int main(int argc, char** argv) { ripser(std::move(dist), dim_max, threshold, ratio, modulus) .compute_barcodes(); } else { - compressed_lower_distance_matrix dist = read_file(filename ? file_stream : std::cin, format); @@ -1057,7 +1032,6 @@ int main(int argc, char** argv) { for (index_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; } @@ -1068,7 +1042,6 @@ int main(int argc, char** argv) { d != std::numeric_limits::infinity() ? std::max(max, d) : max_finite; if (d <= threshold) ++num_edges; } - std::cout << "value range: [" << min << "," << max_finite << "]" << std::endl; if (threshold >= max) { -- cgit v1.2.3 From be5f5cb93477df12b388ff9af7ea1ac403be8356 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Fri, 28 Jun 2019 20:50:56 +0200 Subject: unsigned coeffs --- ripser.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 2455291..9b16200 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -62,7 +62,7 @@ template class hash_map : public std::unordered_map typedef float value_t; typedef int64_t index_t; -typedef int16_t coefficient_t; +typedef uint16_t coefficient_t; class binomial_coeff_table { std::vector> B; @@ -670,7 +670,7 @@ public: for (index_t w : vertices) coface_diameter = std::max(coface_diameter, dist(v, w)); index_t coface_index = idx_above + binomial_coeff(v--, k + 1) + idx_below; coefficient_t coface_coefficient = - (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; + (k & 1 ? modulus - 1 : 1) * get_coefficient(simplex) % modulus; return diameter_entry_t(coface_diameter, coface_index, coface_coefficient); } }; @@ -733,7 +733,7 @@ public: value_t coface_diameter = std::max(get_diameter(simplex), get_diameter(neighbor)); index_t coface_index = idx_above + binomial_coeff(get_index(neighbor), k + 1) + idx_below; coefficient_t coface_coefficient = - (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; + (k & 1 ? modulus - 1 : 1) * get_coefficient(simplex) % modulus; return diameter_entry_t(coface_diameter, coface_index, coface_coefficient); } }; -- cgit v1.2.3 From a75784e9731ef29bde14b15e697f4c112058b46c Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sat, 29 Jun 2019 14:32:38 +0200 Subject: more updates and cleanup --- README.md | 8 +++++--- ripser.cpp | 40 +++++++++++++++++++--------------------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/README.md b/README.md index 3c556c0..46d1ca9 100644 --- a/README.md +++ b/README.md @@ -39,7 +39,7 @@ Ripser's efficiency is based on a few important concepts and principles, buildin ### Version -[Latest release][latest-release]: 1.1 (June 2019) +[Latest release][latest-release]: 1.1 (July 2019) ### Building @@ -89,9 +89,11 @@ The input is given either in a file whose name is passed as an argument, or thro ### Experimental features -The following features are currently planned for future versions: +The following experimental features are currently available in separate branches: - - computation of representative cycles for persistent homology (currenly only *co*cycles are computed) +- `representative-cocycles`: output of representative cocycles for persistent cohomology. +- `representative-cycles`: computation and output of representative cycles for persistent homology (in the standard version, only *co*cycles are computed). +- `simple`: a simplified version of Ripser, without support for sparse distance matrices and coefficients. This might be a good starting point for exploring the code. Prototype implementations are already avaliable; please contact the author if one of these features might be relevant for your research. diff --git a/ripser.cpp b/ripser.cpp index 9b16200..e9ccad7 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -69,13 +69,12 @@ class binomial_coeff_table { public: binomial_coeff_table(index_t n, index_t k) : B(n + 1) { - for (index_t i = 0; i <= n; i++) { - B[i].resize(k + 1); - for (index_t j = 0; j <= std::min(i, k); j++) - if (j == 0 || j == i) - B[i][j] = 1; - else - B[i][j] = B[i - 1][j - 1] + B[i - 1][j]; + for (index_t i = 0; i <= n; ++i) { + B[i].resize(k + 1, 0); + B[i][0] = 1; + if (i <= k) B[i][i] = 1; + for (index_t j = 1; j < std::min(i, k + 1); ++j) + B[i][j] = B[i - 1][j - 1] + B[i - 1][j]; } } @@ -97,7 +96,7 @@ std::vector multiplicative_inverse_vector(const coefficient_t m) inverse[1] = 1; // m = a * (m / a) + m % a // Multipying with inverse(a) * inverse(m % a): - // 0 = inverse(m % a) * (m / a) + inverse(a) (mod m) + // 0 = inverse(m % a) * (m / a) + inverse(a) (mod m) for (coefficient_t a = 2; a < m; ++a) inverse[a] = m - (inverse[m % a] * (m / a)) % m; return inverse; } @@ -395,9 +394,8 @@ public: ratio(_ratio), modulus(_modulus), binomial_coeff(n, dim_max + 2), multiplicative_inverse(multiplicative_inverse_vector(_modulus)) {} - index_t get_next_vertex(index_t& v, const index_t idx, const index_t k) const { - return v = upper_bound( - v, [&](const index_t& w) -> bool { return (binomial_coeff(w, k) <= idx); }); + index_t get_max_vertex(const index_t idx, const index_t k, const index_t n) const { + return upper_bound(n, [&](index_t w) -> bool { return (binomial_coeff(w, k) <= idx); }); } index_t get_edge_index(const index_t i, const index_t j) const { @@ -405,13 +403,13 @@ public: } template - OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v, + OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t n, OutputIterator out) const { - --v; + --n; for (index_t k = dim + 1; k > 0; --k) { - get_next_vertex(v, idx, k); - *out++ = v; - idx -= binomial_coeff(v, k); + n = get_max_vertex(idx, k, n); + *out++ = n; + idx -= binomial_coeff(n, k); } return out; } @@ -501,7 +499,7 @@ public: Column& working_reduction_column, Column& working_coboundary, const index_t& dim, hash_map& pivot_column_index, - bool& might_be_apparent_pair) { + bool& looking_for_emgergent_pair) { 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); @@ -514,10 +512,10 @@ public: diameter_entry_t coface = cofaces.next(); if (get_diameter(coface) <= threshold) { coface_entries.push_back(coface); - if (might_be_apparent_pair && (get_diameter(simplex) == get_diameter(coface))) { + if (looking_for_emgergent_pair && (get_diameter(simplex) == get_diameter(coface))) { if (pivot_column_index.find(get_index(coface)) == pivot_column_index.end()) return coface; - might_be_apparent_pair = false; + looking_for_emgergent_pair = false; } } } @@ -561,14 +559,14 @@ public: std::priority_queue, greater_diameter_or_smaller_index> working_reduction_column, working_coboundary; - bool might_be_apparent_pair = (index_column_to_reduce == index_column_to_add); + bool looking_for_emgergent_pair = (index_column_to_reduce == index_column_to_add); diameter_entry_t pivot; while (true) { pivot = add_coboundary_and_get_pivot(reduction_matrix.cbegin(index_column_to_add), reduction_matrix.cend(index_column_to_add), factor_column_to_add, working_reduction_column, working_coboundary, dim, pivot_column_index, - might_be_apparent_pair); + looking_for_emgergent_pair); if (get_index(pivot) != -1) { auto pair = pivot_column_index.find(get_index(pivot)); if (pair != pivot_column_index.end()) { -- cgit v1.2.3 From aeafba4ba9a3fdfffc84067a262cb0b6bed4a405 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sat, 29 Jun 2019 14:40:36 +0200 Subject: rename check for emergent pair --- ripser.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index e9ccad7..979967e 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -499,7 +499,7 @@ public: Column& working_reduction_column, Column& working_coboundary, const index_t& dim, hash_map& pivot_column_index, - bool& looking_for_emgergent_pair) { + bool& check_for_emergent_pair) { 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); @@ -512,10 +512,11 @@ public: diameter_entry_t coface = cofaces.next(); if (get_diameter(coface) <= threshold) { coface_entries.push_back(coface); - if (looking_for_emgergent_pair && (get_diameter(simplex) == get_diameter(coface))) { + if (check_for_emergent_pair && + (get_diameter(simplex) == get_diameter(coface))) { if (pivot_column_index.find(get_index(coface)) == pivot_column_index.end()) return coface; - looking_for_emgergent_pair = false; + check_for_emergent_pair = false; } } } @@ -559,14 +560,14 @@ public: std::priority_queue, greater_diameter_or_smaller_index> working_reduction_column, working_coboundary; - bool looking_for_emgergent_pair = (index_column_to_reduce == index_column_to_add); + bool check_for_emergent_pair = (index_column_to_reduce == index_column_to_add); diameter_entry_t pivot; while (true) { pivot = add_coboundary_and_get_pivot(reduction_matrix.cbegin(index_column_to_add), reduction_matrix.cend(index_column_to_add), factor_column_to_add, working_reduction_column, working_coboundary, dim, pivot_column_index, - looking_for_emgergent_pair); + check_for_emergent_pair); if (get_index(pivot) != -1) { auto pair = pivot_column_index.find(get_index(pivot)); if (pair != pivot_column_index.end()) { -- cgit v1.2.3 From f7ece0c11c29fe676ea4e03c27035db0fde18e4b Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sat, 29 Jun 2019 23:28:17 +0200 Subject: cleanup and restructuring --- ripser.cpp | 186 ++++++++++++++++++++++++++++++------------------------------- 1 file changed, 93 insertions(+), 93 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 979967e..fde30ad 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -79,7 +79,7 @@ public: } index_t operator()(index_t n, index_t k) const { - assert(n < B.size() && k < B[n].size()); + assert(n < B.size() && k < B[n].size() && n >= k - 1); return B[n][k]; } }; @@ -191,6 +191,35 @@ template struct compressed_distance_matrix { void init_rows(); }; +typedef compressed_distance_matrix compressed_lower_distance_matrix; +typedef compressed_distance_matrix compressed_upper_distance_matrix; + +template <> void compressed_lower_distance_matrix::init_rows() { + value_t* pointer = &distances[0]; + for (index_t i = 1; i < size(); ++i) { + rows[i] = pointer; + pointer += i; + } +} + +template <> void compressed_upper_distance_matrix::init_rows() { + value_t* pointer = &distances[0] - 1; + for (index_t i = 0; i < size() - 1; ++i) { + rows[i] = pointer; + pointer += size() - i - 2; + } +} + +template <> +value_t compressed_lower_distance_matrix::operator()(const index_t i, const index_t j) const { + return i == j ? 0 : i < j ? rows[j][i] : rows[i][j]; +} + +template <> +value_t compressed_upper_distance_matrix::operator()(const index_t i, const index_t j) const { + return i == j ? 0 : i > j ? rows[j][i] : rows[i][j]; +} + struct sparse_distance_matrix { std::vector> neighbors; @@ -218,37 +247,6 @@ struct sparse_distance_matrix { size_t size() const { return neighbors.size(); } }; -template <> void compressed_distance_matrix::init_rows() { - value_t* pointer = &distances[0]; - for (index_t i = 1; i < size(); ++i) { - rows[i] = pointer; - pointer += i; - } -} - -template <> void compressed_distance_matrix::init_rows() { - value_t* pointer = &distances[0] - 1; - for (index_t i = 0; i < size() - 1; ++i) { - rows[i] = pointer; - pointer += size() - i - 2; - } -} - -template <> -value_t compressed_distance_matrix::operator()(const index_t i, - const index_t j) const { - return i == j ? 0 : i > j ? rows[j][i] : rows[i][j]; -} - -template <> -value_t compressed_distance_matrix::operator()(const index_t i, - const index_t j) const { - return i == j ? 0 : i < j ? rows[j][i] : rows[i][j]; -} - -typedef compressed_distance_matrix compressed_lower_distance_matrix; -typedef compressed_distance_matrix compressed_upper_distance_matrix; - struct euclidean_distance_matrix { std::vector> points; @@ -298,7 +296,7 @@ public: } }; -template diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) { +template diameter_entry_t pop_pivot(Column& column, coefficient_t modulus) { diameter_entry_t pivot(-1); while (!column.empty()) { if (get_coefficient(pivot) == 0) @@ -314,7 +312,7 @@ template diameter_entry_t pop_pivot(Heap& column, coefficient_t return pivot; } -template diameter_entry_t get_pivot(Heap& column, coefficient_t modulus) { +template diameter_entry_t get_pivot(Column& column, coefficient_t modulus) { diameter_entry_t result = pop_pivot(column, modulus); if (get_index(result) != -1) column.push(result); return result; @@ -327,7 +325,7 @@ template class compressed_sparse_matrix { public: size_t size() const { return bounds.size(); } - typename std::vector::const_iterator cbegin(size_t index) const { + typename std::vector::const_iterator cbegin(const size_t index) const { assert(index < size()); return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1]; } @@ -337,7 +335,7 @@ public: return entries.cbegin() + bounds[index]; } - template void append_column(Iterator begin, Iterator end) { + template void append_column(const Iterator begin, const Iterator end) { for (Iterator it = begin; it != end; ++it) { entries.push_back(*it); } bounds.push_back(entries.size()); } @@ -361,13 +359,13 @@ public: } }; -template index_t upper_bound(index_t top, Predicate pred) { +template index_t get_max(index_t top, index_t bottom, const Predicate pred) { if (!pred(top)) { - index_t count = top; + index_t count = top - bottom; while (count > 0) { - index_t step = count >> 1; - if (!pred(top - step)) { - top -= step + 1; + index_t step = count >> 1, mid = top - step; + if (!pred(mid)) { + top = mid - 1; count -= step + 1; } else count = step; @@ -377,13 +375,13 @@ template index_t upper_bound(index_t top, Predicate pred) { } template class ripser { - DistanceMatrix dist; - index_t n, dim_max; - value_t threshold; - float ratio; - coefficient_t modulus; + const DistanceMatrix dist; + const index_t n, dim_max; + const value_t threshold; + const float ratio; + const coefficient_t modulus; const binomial_coeff_table binomial_coeff; - std::vector multiplicative_inverse; + const std::vector multiplicative_inverse; mutable std::vector coface_entries; public: @@ -395,7 +393,7 @@ public: multiplicative_inverse(multiplicative_inverse_vector(_modulus)) {} index_t get_max_vertex(const index_t idx, const index_t k, const index_t n) const { - return upper_bound(n, [&](index_t w) -> bool { return (binomial_coeff(w, k) <= idx); }); + return get_max(n, k - 1, [&](index_t w) -> bool { return (binomial_coeff(w, k) <= idx); }); } index_t get_edge_index(const index_t i, const index_t j) const { @@ -493,36 +491,44 @@ public: #endif } + 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) { + bool check_for_emergent_pair = true; + coface_entries.clear(); + 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 (check_for_emergent_pair && (get_diameter(simplex) == get_diameter(coface))) { + if (pivot_column_index.find(get_index(coface)) == pivot_column_index.end()) + return coface; + check_for_emergent_pair = false; + } + } + } + for (auto coface : coface_entries) working_coboundary.push(coface); + return get_pivot(working_coboundary, modulus); + } + template - diameter_entry_t add_coboundary_and_get_pivot(Iterator column_begin, Iterator column_end, - coefficient_t factor_column_to_add, - Column& working_reduction_column, - Column& working_coboundary, const index_t& dim, - hash_map& pivot_column_index, - bool& check_for_emergent_pair) { + 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(); 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); - working_reduction_column.push(simplex); - - coface_entries.clear(); 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 (check_for_emergent_pair && - (get_diameter(simplex) == get_diameter(coface))) { - if (pivot_column_index.find(get_index(coface)) == pivot_column_index.end()) - return coface; - check_for_emergent_pair = false; - } - } + if (get_diameter(coface) <= threshold) coface_entries.push_back(coface); } - for (auto coface : coface_entries) working_coboundary.push(coface); } - return get_pivot(working_coboundary, modulus); + for (auto coface : coface_entries) working_coboundary.push(coface); } void compute_pairs(std::vector& columns_to_reduce, @@ -536,7 +542,8 @@ public: for (index_t index_column_to_reduce = 0; index_column_to_reduce < columns_to_reduce.size(); ++index_column_to_reduce) { - auto column_to_reduce = columns_to_reduce[index_column_to_reduce]; + + diameter_entry_t column_to_reduce(columns_to_reduce[index_column_to_reduce], 1); value_t diameter = get_diameter(column_to_reduce); #ifdef INDICATE_PROGRESS @@ -547,32 +554,28 @@ public: << "\r" << std::flush; #endif - index_t index_column_to_add = index_column_to_reduce; - - // start with factor 1 in order to initialize working_coboundary - // with the coboundary of the simplex with index column_to_reduce - coefficient_t factor_column_to_add = 1; - - // initialize reduction_matrix as identity matrix - reduction_matrix.append_column(); - reduction_matrix.push_back(diameter_entry_t(column_to_reduce, 1)); - std::priority_queue, greater_diameter_or_smaller_index> working_reduction_column, working_coboundary; - bool check_for_emergent_pair = (index_column_to_reduce == index_column_to_add); - diameter_entry_t pivot; + + 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); + while (true) { - pivot = add_coboundary_and_get_pivot(reduction_matrix.cbegin(index_column_to_add), - reduction_matrix.cend(index_column_to_add), - factor_column_to_add, working_reduction_column, - working_coboundary, dim, pivot_column_index, - check_for_emergent_pair); if (get_index(pivot) != -1) { auto pair = pivot_column_index.find(get_index(pivot)); if (pair != pivot_column_index.end()) { - index_column_to_add = pair->second; - factor_column_to_add = modulus - get_coefficient(pivot); + index_t index_column_to_add = pair->second; + coefficient_t factor_column_to_add = modulus - get_coefficient(pivot); + + 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); + + pivot = get_pivot(working_coboundary, modulus); } else { #ifdef PRINT_PERSISTENCE_PAIRS value_t death = get_diameter(pivot); @@ -586,10 +589,7 @@ public: pivot_column_index.insert( std::make_pair(get_index(pivot), index_column_to_reduce)); - // replace current column of reduction_matrix (with a single diagonal 1 - // entry) by reduction_column (possibly with a different entry on the - // diagonal) - reduction_matrix.pop_back(); + reduction_matrix.append_column(); const coefficient_t inverse = multiplicative_inverse[get_coefficient(pivot)]; @@ -644,7 +644,7 @@ template <> class ripser::simplex_coboundary_e const binomial_coeff_table& binomial_coeff; 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) : 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), -- cgit v1.2.3 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 From dcec6219f9bb245fa265063c6dd1c9d06d15aab0 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Tue, 2 Jul 2019 11:51:43 +0200 Subject: fixed misplaced call to reduction_matrix.append_column --- ripser.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index fde30ad..02527d4 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -554,6 +554,8 @@ public: << "\r" << std::flush; #endif + reduction_matrix.append_column(); + std::priority_queue, greater_diameter_or_smaller_index> working_reduction_column, working_coboundary; @@ -589,8 +591,6 @@ public: pivot_column_index.insert( std::make_pair(get_index(pivot), index_column_to_reduce)); - reduction_matrix.append_column(); - const coefficient_t inverse = multiplicative_inverse[get_coefficient(pivot)]; while (true) { -- cgit v1.2.3 From b909e6c124b8d83e759793991f9db391103ebb0b Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Tue, 2 Jul 2019 12:52:59 +0200 Subject: cleanup add_coboundary --- benchmarks/Dockerfile | 16 +++++++++++----- ripser.cpp | 21 +++++++++++---------- 2 files changed, 22 insertions(+), 15 deletions(-) diff --git a/benchmarks/Dockerfile b/benchmarks/Dockerfile index 40c654f..3d6e952 100644 --- a/benchmarks/Dockerfile +++ b/benchmarks/Dockerfile @@ -52,7 +52,7 @@ RUN curl -LO https://raw.githubusercontent.com/Ripser/ripser/dev/benchmarks/sphe && curl -LO https://github.com/n-otter/PH-roadmap/raw/master/data_sets/roadmap_datasets_point_cloud/random_point_cloud_50_16_.txt \ && curl -LO https://github.com/n-otter/PH-roadmap/raw/master/data_sets/roadmap_datasets_distmat/dipha/random_50_16.bin - +COPY o3_1024.txt . FROM benchmark-setup as benchmark-ripser @@ -65,10 +65,11 @@ ENV PATH="${PATH}:/ripser/ripser" WORKDIR /benchmark RUN time -v -o sphere_3_192.ripser.txt ripser sphere_3_192.distance_matrix --dim 2 -RUN time -v -o random.ripser.txt ripser random_point_cloud_50_16_.txt --dim 7 +RUN time -v -o random.ripser.txt ripser random_point_cloud_50_16_.txt --format point-cloud --dim 7 RUN time -v -o fractal-r.ripser.txt ripser fractal_9_5_2_random_edge_list.txt_0.19795_distmat.txt --dim 2 -RUN time -v -o dragon-2.ripser.txt ripser dragon_vrip.ply.txt_2000_.txt --dim 1 +RUN time -v -o dragon-2.ripser.txt ripser dragon_vrip.ply.txt_2000_.txt --format point-cloud --dim 1 # RUN time -v -o genome.ripser.txt ripser human_gene_distmat.txt --dim 1 +RUN time -v -o o3_1024.ripser.txt ripser o3_1024.txt --format point-cloud --dim 3 --threshold 1.8 --ratio 2 @@ -175,8 +176,13 @@ RUN time -v -o fractal-r.eirene.txt \ julia --eval 'using Eirene; using DelimitedFiles; eirene([0 1; 1 0]); println(@elapsed eirene(transpose(readdlm("fractal_9_5_2_random_edge_list.txt_0.19795_distmat.txt")), record="none", maxdim=2))' >fractal-r.eirene.out.txt #RUN time -v -o random.eirene.txt \ # julia --eval 'using Eirene; using DelimitedFiles; eirene([0 1; 1 0]); println(@elapsed eirene(transpose(readdlm("random_point_cloud_50_16_.txt")), model="pc", record="none", maxdim=7))' >random.eirene.out.txt - - +RUN time -v -o o3_1024.eirene.txt \ + julia --eval 'using Eirene; using DelimitedFiles; eirene([0 1; 1 0]); println(@elapsed eirene(transpose(readdlm("o3_1024.txt")), model="pc", record="none", maxdim=3, maxrad=1.8))' > o3_1024.eirene.out.txt +# +# using Eirene; using DelimitedFiles; eirene([0 1; 1 0]); println(@elapsed result=eirene(transpose(readdlm("/Users/uli/Bitbucket/ripser/examples/o3_1024.txt")), model="pc", record="none", maxdim=3, maxrad=1.8)) +# using Eirene; using DelimitedFiles; eirene([0 1; 1 0]); println(@elapsed result=eirene(transpose(readdlm("/Users/uli/Bitbucket/ripser/examples/o3_2048.txt")), model="pc", record="none", maxdim=3, maxrad=1.6)) +# using Eirene; using DelimitedFiles; eirene([0 1; 1 0]); println(@elapsed result=eirene(transpose(readdlm("/Users/uli/Bitbucket/ripser/examples/o3_4096.txt")), model="pc", record="none", maxdim=3, maxrad=1.4)) +# FROM benchmark-setup as benchmark-dionysus2 diff --git a/ripser.cpp b/ripser.cpp index 3f19858..b10ec07 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -539,12 +539,18 @@ public: if (get_diameter(coface) <= threshold) working_coboundary.push(coface); } } - + template void add_coboundary(compressed_sparse_matrix& reduction_matrix, - const index_t index_column_to_add, const coefficient_t factor, - Column& working_reduction_column, Column& working_coboundary, - const index_t& dim) { + std::vector& columns_to_reduce, + const index_t index_column_to_add, const coefficient_t factor, + Column& working_reduction_column, Column& working_coboundary, + const index_t& dim) { + 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); + for (auto it = reduction_matrix.cbegin(index_column_to_add); it != reduction_matrix.cend(index_column_to_add); ++it) { diameter_entry_t simplex = *it; @@ -598,12 +604,7 @@ 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); - - add_coboundary(column_to_add, working_reduction_column, working_coboundary, - dim); - add_coboundary(reduction_matrix, index_column_to_add, factor, + add_coboundary(reduction_matrix, columns_to_reduce, index_column_to_add, factor, working_reduction_column, working_coboundary, dim); pivot = get_pivot(working_coboundary, modulus); -- cgit v1.2.3 From 74e20affd9256b78467729e9f9d9810abbd36c9d Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Tue, 2 Jul 2019 20:43:26 +0200 Subject: changed interface to compressed_sparse_matrix --- ripser.cpp | 45 ++++++++++++++++++++++----------------------- 1 file changed, 22 insertions(+), 23 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index b10ec07..0a32a6d 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -333,6 +333,13 @@ template diameter_entry_t get_pivot(Column& column, const coef return result; } +template struct iterator_range { + Iterator b, e; + iterator_range(Iterator _b, Iterator _e) : b(_b), e(_e) {} + Iterator begin() const { return b; } + Iterator end() const { return e; } +}; + template class compressed_sparse_matrix { std::vector bounds; std::vector entries; @@ -345,14 +352,10 @@ public: return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1]; } - typename std::vector::const_iterator cend(const size_t index) const { - assert(index < size()); - return entries.cbegin() + bounds[index]; - } - - template void append_column(const Iterator begin, const Iterator end) { - for (Iterator it = begin; it != end; ++it) { entries.push_back(*it); } - bounds.push_back(entries.size()); + iterator_range::const_iterator> subrange(const index_t index) { + return iterator_range::const_iterator>( + index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1], + entries.cbegin() + bounds[index]); } void append_column() { bounds.push_back(entries.size()); } @@ -539,21 +542,17 @@ public: if (get_diameter(coface) <= threshold) working_coboundary.push(coface); } } - + template void add_coboundary(compressed_sparse_matrix& reduction_matrix, - std::vector& columns_to_reduce, - const index_t index_column_to_add, const coefficient_t factor, - Column& working_reduction_column, Column& working_coboundary, - const index_t& dim) { - 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); - - for (auto it = reduction_matrix.cbegin(index_column_to_add); - it != reduction_matrix.cend(index_column_to_add); ++it) { - diameter_entry_t simplex = *it; + std::vector& columns_to_reduce, + const index_t index_column_to_add, const coefficient_t factor, + Column& working_reduction_column, Column& working_coboundary, + const index_t& dim) { + 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); + + for (auto simplex : reduction_matrix.subrange(index_column_to_add)) { set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); working_reduction_column.push(simplex); add_coboundary(simplex, working_reduction_column, working_coboundary, dim); @@ -604,8 +603,8 @@ public: modulus - get_coefficient(pivot) * multiplicative_inverse[get_coefficient(other_pivot)] % modulus; - add_coboundary(reduction_matrix, columns_to_reduce, index_column_to_add, factor, - working_reduction_column, working_coboundary, dim); + add_coboundary(reduction_matrix, columns_to_reduce, index_column_to_add, + factor, working_reduction_column, working_coboundary, dim); pivot = get_pivot(working_coboundary, modulus); } else { -- cgit v1.2.3 From 028b237e48d92efd9655a920ed38d0c441540ce8 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Tue, 2 Jul 2019 20:45:27 +0200 Subject: some cleanup --- ripser.cpp | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index b10ec07..3d1a9e9 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -539,18 +539,16 @@ public: if (get_diameter(coface) <= threshold) working_coboundary.push(coface); } } - + template void add_coboundary(compressed_sparse_matrix& reduction_matrix, - std::vector& columns_to_reduce, - const index_t index_column_to_add, const coefficient_t factor, - Column& working_reduction_column, Column& working_coboundary, - const index_t& dim) { - 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); - + const std::vector& columns_to_reduce, + const index_t index_column_to_add, const coefficient_t factor, + Column& working_reduction_column, Column& working_coboundary, + const index_t& dim) { + 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); + for (auto it = reduction_matrix.cbegin(index_column_to_add); it != reduction_matrix.cend(index_column_to_add); ++it) { diameter_entry_t simplex = *it; @@ -560,7 +558,7 @@ public: } } - void compute_pairs(std::vector& columns_to_reduce, + void compute_pairs(const std::vector& columns_to_reduce, entry_hash_map& pivot_column_index, const index_t dim) { #ifdef PRINT_PERSISTENCE_PAIRS @@ -604,8 +602,9 @@ public: modulus - get_coefficient(pivot) * multiplicative_inverse[get_coefficient(other_pivot)] % modulus; - add_coboundary(reduction_matrix, columns_to_reduce, index_column_to_add, factor, - working_reduction_column, working_coboundary, dim); + + add_coboundary(reduction_matrix, columns_to_reduce, index_column_to_add, + factor, working_reduction_column, working_coboundary, dim); pivot = get_pivot(working_coboundary, modulus); } else { -- cgit v1.2.3 From 36953ea242ef98a14809978bbe2f04aa96f981ca Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 3 Jul 2019 16:29:12 +0200 Subject: cleanup --- ripser.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 0822eb2..422e850 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -588,8 +588,6 @@ public: greater_diameter_or_smaller_index> working_reduction_column, working_coboundary; - // 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); @@ -621,8 +619,6 @@ public: pivot_column_index.insert( std::make_pair(get_entry(pivot), index_column_to_reduce)); - 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; -- cgit v1.2.3 From 6c6a6ef9dbecc3bca4c67db5f8dd9af8eb68e1b9 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 3 Jul 2019 21:24:54 +0200 Subject: cleanup --- ripser.cpp | 54 +++++++++++++++++------------------------------------- 1 file changed, 17 insertions(+), 37 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 422e850..08b3f2a 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -114,9 +114,6 @@ struct __attribute__((packed)) entry_t { static_assert(sizeof(entry_t) == sizeof(index_t), "size of entry_t is not the same as index_t"); -entry_t make_entry(index_t _index, coefficient_t _coefficient) { - return entry_t(_index, _coefficient); -} index_t get_index(const entry_t& e) { return e.index; } index_t get_coefficient(const entry_t& e) { return e.coefficient; } void set_coefficient(entry_t& e, const coefficient_t c) { e.coefficient = c; } @@ -153,12 +150,12 @@ index_t get_index(const index_diameter_t& i) { return i.first; } value_t get_diameter(const index_diameter_t& i) { return i.second; } struct diameter_entry_t : std::pair { - diameter_entry_t() {} + using std::pair::pair; diameter_entry_t(value_t _diameter, index_t _index, coefficient_t _coefficient) - : std::pair(_diameter, make_entry(_index, _coefficient)) {} + : diameter_entry_t(_diameter, {_index, _coefficient}) {} diameter_entry_t(const diameter_index_t& _diameter_index, coefficient_t _coefficient) - : std::pair(get_diameter(_diameter_index), - make_entry(get_index(_diameter_index), _coefficient)) {} + : diameter_entry_t(get_diameter(_diameter_index), + {get_index(_diameter_index), _coefficient}) {} diameter_entry_t(const index_t& _index) : diameter_entry_t(0, _index, 0) {} }; @@ -255,7 +252,7 @@ struct sparse_distance_matrix { for (index_t j = 0; j < size(); ++j) if (i != j && mat(i, j) <= threshold) { ++num_edges; - neighbors[i].push_back(std::make_pair(j, mat(i, j))); + neighbors[i].push_back({j, mat(i, j)}); } } @@ -333,29 +330,22 @@ template diameter_entry_t get_pivot(Column& column, const coef return result; } -template struct iterator_range { - Iterator b, e; - iterator_range(Iterator _b, Iterator _e) : b(_b), e(_e) {} - Iterator begin() const { return b; } - Iterator end() const { return e; } -}; +template T begin(std::pair& p) { return p.first; } +template T end(std::pair& p) { return p.second; } template class compressed_sparse_matrix { std::vector bounds; std::vector entries; + typedef typename std::vector::iterator iterator; + typedef std::pair iterator_pair; + public: size_t size() const { return bounds.size(); } - typename std::vector::const_iterator cbegin(const size_t index) const { - assert(index < size()); - return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1]; - } - - iterator_range::const_iterator> subrange(const index_t index) { - return iterator_range::const_iterator>( - index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1], - entries.cbegin() + bounds[index]); + iterator_pair subrange(const index_t index) { + return {entries.begin() + (index == 0 ? 0 : bounds[index - 1]), + entries.begin() + bounds[index]}; } void append_column() { bounds.push_back(entries.size()); } @@ -365,16 +355,6 @@ public: entries.push_back(e); ++bounds.back(); } - - void pop_back() { - assert(0 < size()); - entries.pop_back(); - --bounds.back(); - } - - template void append_column(const Collection collection) { - append_column(collection.cbegin(), collection.cend()); - } }; template @@ -768,7 +748,7 @@ template <> std::vector ripser 0;) { get_simplex_vertices(index, 1, dist.size(), vertices.rbegin()); value_t length = dist(vertices[0], vertices[1]); - if (length <= threshold) edges.push_back(std::make_pair(length, index)); + if (length <= threshold) edges.push_back({length, index}); } return edges; } @@ -778,7 +758,7 @@ template <> std::vector ripser::get_ed for (index_t i = 0; i < n; ++i) for (auto n : dist.neighbors[i]) { index_t j = get_index(n); - if (i > j) edges.push_back(std::make_pair(get_diameter(n), get_edge_index(i, j))); + if (i > j) edges.push_back({get_diameter(n), get_edge_index(i, j)}); } return edges; } @@ -841,8 +821,8 @@ sparse_distance_matrix read_sparse_distance_matrix(std::istream& input_stream) { s >> value; if (i != j) { neighbors.resize(std::max({neighbors.size(), i + 1, j + 1})); - neighbors[i].push_back(std::make_pair(j, value)); - neighbors[j].push_back(std::make_pair(i, value)); + neighbors[i].push_back({j, value}); + neighbors[j].push_back({i, value}); ++num_edges; } } -- cgit v1.2.3 From 7bbb01b86defc3314aa32ed5284104039a82fc54 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 3 Jul 2019 21:32:27 +0200 Subject: better progress indication --- ripser.cpp | 62 ++++++++++++++++++++++++++++++++++++++++---------------------- 1 file changed, 40 insertions(+), 22 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 08b3f2a..75af26a 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -36,7 +36,7 @@ */ -//#define INDICATE_PROGRESS +#define INDICATE_PROGRESS #define PRINT_PERSISTENCE_PAIRS //#define USE_GOOGLE_HASHMAP @@ -47,6 +47,7 @@ #include #include #include +#include #ifdef USE_GOOGLE_HASHMAP #include @@ -66,6 +67,8 @@ typedef float value_t; typedef int64_t index_t; typedef uint16_t coefficient_t; +static const std::chrono::milliseconds time_step(40); + class binomial_coeff_table { std::vector> B; @@ -418,42 +421,55 @@ public: entry_hash_map& pivot_column_index, index_t dim) { #ifdef INDICATE_PROGRESS - std::cout << "\033[K" + std::cerr << "\r\033[K" << "assembling columns" - << "\r" << std::flush; + << std::flush; #endif + std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step; --dim; columns_to_reduce.clear(); std::vector next_simplices; + + index_t i = 0; - for (diameter_index_t simplex : simplices) { + for (diameter_index_t& simplex : simplices) { + ++i; + simplex_coboundary_enumerator cofaces(diameter_entry_t(simplex, 1), dim, *this); while (cofaces.has_next(false)) { +#ifdef INDICATE_PROGRESS + if (std::chrono::steady_clock::now() > next) { + std::cerr << "\r\033[K" + << "assembling " << next_simplices.size() << " columns (processing " + << std::distance(&simplices[0], &simplex) << "/" << simplices.size() << " simplices)" + << std::flush; + next = std::chrono::steady_clock::now() + time_step; + } +#endif auto coface = cofaces.next(); - next_simplices.push_back(std::make_pair(get_diameter(coface), get_index(coface))); + next_simplices.push_back({get_diameter(coface), get_index(coface)}); 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))); + columns_to_reduce.push_back({get_diameter(coface), get_index(coface)}); } } simplices.swap(next_simplices); #ifdef INDICATE_PROGRESS - std::cout << "\033[K" + std::cerr << "\r\033[K" << "sorting " << columns_to_reduce.size() << " columns" - << "\r" << std::flush; + << std::flush; #endif std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), greater_diameter_or_smaller_index()); #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cerr << "\r\033[K"; #endif } @@ -547,6 +563,8 @@ public: #endif compressed_sparse_matrix reduction_matrix; + + std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step; for (index_t index_column_to_reduce = 0; index_column_to_reduce < columns_to_reduce.size(); ++index_column_to_reduce) { @@ -554,14 +572,6 @@ public: diameter_entry_t column_to_reduce(columns_to_reduce[index_column_to_reduce], 1); value_t diameter = get_diameter(column_to_reduce); -#ifdef INDICATE_PROGRESS - if ((index_column_to_reduce + 1) % 1000000 == 0) - std::cout << "\033[K" - << "reducing column " << index_column_to_reduce + 1 << "/" - << columns_to_reduce.size() << " (diameter " << diameter << ")" - << "\r" << std::flush; -#endif - reduction_matrix.append_column(); std::priority_queue, @@ -572,6 +582,15 @@ public: column_to_reduce, working_coboundary, dim, pivot_column_index); while (true) { +#ifdef INDICATE_PROGRESS + if (std::chrono::steady_clock::now() > next) { + std::cerr << "\r\033[K" + << "reducing column " << index_column_to_reduce + 1 << "/" + << columns_to_reduce.size() << " (diameter " << diameter << ")" + << std::flush; + next = std::chrono::steady_clock::now() + time_step; + } +#endif if (get_index(pivot) != -1) { auto pair = pivot_column_index.find(get_entry(pivot)); if (pair != pivot_column_index.end()) { @@ -591,13 +610,12 @@ public: value_t death = get_diameter(pivot); if (death > diameter * ratio) { #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cerr << "\r\033[K"; #endif 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({get_entry(pivot), index_column_to_reduce}); while (true) { diameter_entry_t e = pop_pivot(working_reduction_column, modulus); @@ -616,7 +634,7 @@ public: } } #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cerr << "\r\033[K" << std::flush; #endif } -- cgit v1.2.3 From e843303d304ffff9b1fa8e8c073ceaa50469fe59 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 3 Jul 2019 23:24:12 +0200 Subject: clear line --- ripser.cpp | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 75af26a..eedfa2b 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -69,6 +69,8 @@ typedef uint16_t coefficient_t; static const std::chrono::milliseconds time_step(40); +static const std::string clear_line("\r\033[K"); + class binomial_coeff_table { std::vector> B; @@ -421,7 +423,7 @@ public: entry_hash_map& pivot_column_index, index_t dim) { #ifdef INDICATE_PROGRESS - std::cerr << "\r\033[K" + std::cerr << clear_line << "assembling columns" << std::flush; #endif @@ -442,7 +444,7 @@ public: while (cofaces.has_next(false)) { #ifdef INDICATE_PROGRESS if (std::chrono::steady_clock::now() > next) { - std::cerr << "\r\033[K" + std::cerr << clear_line << "assembling " << next_simplices.size() << " columns (processing " << std::distance(&simplices[0], &simplex) << "/" << simplices.size() << " simplices)" << std::flush; @@ -461,7 +463,7 @@ public: simplices.swap(next_simplices); #ifdef INDICATE_PROGRESS - std::cerr << "\r\033[K" + std::cerr << clear_line << "sorting " << columns_to_reduce.size() << " columns" << std::flush; #endif @@ -469,7 +471,7 @@ public: std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), greater_diameter_or_smaller_index()); #ifdef INDICATE_PROGRESS - std::cerr << "\r\033[K"; + std::cerr << clear_line << std::flush; #endif } @@ -584,7 +586,7 @@ public: while (true) { #ifdef INDICATE_PROGRESS if (std::chrono::steady_clock::now() > next) { - std::cerr << "\r\033[K" + std::cerr << clear_line << "reducing column " << index_column_to_reduce + 1 << "/" << columns_to_reduce.size() << " (diameter " << diameter << ")" << std::flush; @@ -610,7 +612,7 @@ public: value_t death = get_diameter(pivot); if (death > diameter * ratio) { #ifdef INDICATE_PROGRESS - std::cerr << "\r\033[K"; + std::cerr << clear_line << std::flush; #endif std::cout << " [" << diameter << "," << death << ")" << std::endl; } @@ -627,6 +629,9 @@ public: } } else { #ifdef PRINT_PERSISTENCE_PAIRS +#ifdef INDICATE_PROGRESS + std::cerr << clear_line << std::flush; +#endif std::cout << " [" << diameter << ", )" << std::endl; #endif break; @@ -634,7 +639,7 @@ public: } } #ifdef INDICATE_PROGRESS - std::cerr << "\r\033[K" << std::flush; + std::cerr << clear_line << std::flush; #endif } -- cgit v1.2.3