From 4e0a06ed43493331ef0e2f5cfe86ed47a7b514ee Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 11 Aug 2019 16:25:31 +0200 Subject: renamed variables, changing "face" to "facet" --- ripser.cpp | 70 +++++++++++++++++++++++++++++++------------------------------- 1 file changed, 35 insertions(+), 35 deletions(-) (limited to 'ripser.cpp') diff --git a/ripser.cpp b/ripser.cpp index 8c345cc..8fb0cd1 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -386,7 +386,7 @@ template class ripser { const coefficient_t modulus; const binomial_coeff_table binomial_coeff; const std::vector multiplicative_inverse; - mutable std::vector coface_entries; + mutable std::vector cofacet_entries; struct entry_hash { std::size_t operator()(const entry_t& e) const { @@ -446,9 +446,9 @@ public: std::vector next_simplices; for (diameter_index_t& simplex : simplices) { - simplex_coboundary_enumerator cofaces(diameter_entry_t(simplex, 1), dim, *this); + simplex_coboundary_enumerator cofacets(diameter_entry_t(simplex, 1), dim, *this); - while (cofaces.has_next(false)) { + while (cofacets.has_next(false)) { #ifdef INDICATE_PROGRESS if (std::chrono::steady_clock::now() > next) { std::cerr << clear_line << "assembling " << next_simplices.size() @@ -457,13 +457,13 @@ public: next = std::chrono::steady_clock::now() + time_step; } #endif - auto coface = cofaces.next(); - if (get_diameter(coface) <= threshold) { + auto cofacet = cofacets.next(); + if (get_diameter(cofacet) <= threshold) { - next_simplices.push_back({get_diameter(coface), get_index(coface)}); + next_simplices.push_back({get_diameter(cofacet), get_index(cofacet)}); - if (pivot_column_index.find(get_entry(coface)) == pivot_column_index.end()) - columns_to_reduce.push_back({get_diameter(coface), get_index(coface)}); + if (pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()) + columns_to_reduce.push_back({get_diameter(cofacet), get_index(cofacet)}); } } } @@ -552,20 +552,20 @@ public: 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); - 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_entry(coface)) == pivot_column_index.end()) - return coface; + cofacet_entries.clear(); + simplex_coboundary_enumerator cofacets(simplex, dim, *this); + while (cofacets.has_next()) { + diameter_entry_t cofacet = cofacets.next(); + if (get_diameter(cofacet) <= threshold) { + cofacet_entries.push_back(cofacet); + if (check_for_emergent_pair && (get_diameter(simplex) == get_diameter(cofacet))) { + if (pivot_column_index.find(get_entry(cofacet)) == pivot_column_index.end()) + return cofacet; check_for_emergent_pair = false; } } } - for (auto coface : coface_entries) working_coboundary.push(coface); + for (auto cofacet : cofacet_entries) working_coboundary.push(cofacet); return get_pivot(working_coboundary); } @@ -573,10 +573,10 @@ public: 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()) { - diameter_entry_t coface = cofaces.next(); - if (get_diameter(coface) <= threshold) working_coboundary.push(coface); + simplex_coboundary_enumerator cofacets(simplex, dim, *this); + while (cofacets.has_next()) { + diameter_entry_t cofacet = cofacets.next(); + if (get_diameter(cofacet) <= threshold) working_coboundary.push(cofacet); } } @@ -719,9 +719,9 @@ public: parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.rbegin()); } - bool has_next(bool all_cofaces = true) { + bool has_next(bool all_cofacets = true) { while ((v != -1) && (binomial_coeff(v, k) <= idx_below)) { - if (!all_cofaces) return false; + if (!all_cofacets) return false; idx_below -= binomial_coeff(v, k); idx_above += binomial_coeff(v, k + 1); --v; @@ -732,12 +732,12 @@ public: } diameter_entry_t next() { - value_t coface_diameter = get_diameter(simplex); - 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 = + value_t cofacet_diameter = get_diameter(simplex); + for (index_t w : vertices) cofacet_diameter = std::max(cofacet_diameter, dist(v, w)); + index_t cofacet_index = idx_above + binomial_coeff(v--, k + 1) + idx_below; + coefficient_t cofacet_coefficient = (k & 1 ? modulus - 1 : 1) * get_coefficient(simplex) % modulus; - return diameter_entry_t(coface_diameter, coface_index, coface_coefficient); + return diameter_entry_t(cofacet_diameter, cofacet_index, cofacet_coefficient); } }; @@ -770,7 +770,7 @@ public: } } - bool has_next(bool all_cofaces = true) { + bool has_next(bool all_cofacets = true) { for (auto &it0 = neighbor_it[0], &end0 = neighbor_end[0]; it0 != end0; ++it0) { neighbor = *it0; for (size_t idx = 1; idx < neighbor_it.size(); ++idx) { @@ -783,7 +783,7 @@ public: neighbor = std::max(neighbor, *it); } while (k > 0 && vertices[k - 1] > get_index(neighbor)) { - if (!all_cofaces) return false; + if (!all_cofacets) return false; idx_below -= binomial_coeff(vertices[k - 1], k); idx_above += binomial_coeff(vertices[k - 1], k + 1); --k; @@ -796,11 +796,11 @@ public: diameter_entry_t next() { ++neighbor_it[0]; - 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 = + value_t cofacet_diameter = std::max(get_diameter(simplex), get_diameter(neighbor)); + index_t cofacet_index = idx_above + binomial_coeff(get_index(neighbor), k + 1) + idx_below; + coefficient_t cofacet_coefficient = (k & 1 ? modulus - 1 : 1) * get_coefficient(simplex) % modulus; - return diameter_entry_t(coface_diameter, coface_index, coface_coefficient); + return diameter_entry_t(cofacet_diameter, cofacet_index, cofacet_coefficient); } }; -- cgit v1.2.3 From a304cb663b2e7e7680c15f8acf5aa36ad9bb9f15 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 11 Aug 2019 23:30:31 +0200 Subject: updated and cleaned up readme --- README.md | 25 ++++++++++++------------- ripser.cpp | 8 ++++---- 2 files changed, 16 insertions(+), 17 deletions(-) (limited to 'ripser.cpp') diff --git a/README.md b/README.md index 445d8eb..36b041c 100644 --- a/README.md +++ b/README.md @@ -24,22 +24,22 @@ 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 + - 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: +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: - Compute persistent *co*homology (as suggested by [Vin de Silva, Dmitriy Morozov, and Mikael Vejdemo-Johansson](https://doi.org/10.1088/0266-5611/27/12/124003)) - - Don't compute information that is never needed - (for the experts: employ the *clearing* optimization, aka *persistence with a twist*, as suggested by [Chao Chen and Michael Kerber](http://www.geometrie.tugraz.at/kerber/kerber_papers/ck-phcwat-11.pdf)) - - Don't store information that can be readily recomputed (in particular, the boundary matrix and the reduced boundary matrix) - - Take computational shortcuts (*apparent* and *emergent persistence pairs*) + - Use the chain complex property that boundaries are cycles + (employ the *clearing* optimization, aka *persistence with a twist*, as suggested by [Chao Chen and Michael Kerber](http://www.geometrie.tugraz.at/kerber/kerber_papers/ck-phcwat-11.pdf)) - If no threshold is specified, choose the *enclosing radius* as the threshold, from which on homology is guaranteed to be trivial (as suggested by [Greg Henselman-Petrusek](https://github.com/Eetion/Eirene.jl)) + - Don't store information that can be readily recomputed (in particular, the original and the reduced boundary matrix) + - Take computational shortcuts (*apparent* and *emergent persistence pairs*) ### Version -[Latest release][latest-release]: 1.1 (July 2019) +[Latest release][latest-release]: 1.1 (August 2019) ### Building @@ -58,6 +58,7 @@ 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: + - `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 @@ -77,12 +78,12 @@ The input is given either in a file whose name is passed as an argument, or thro - `upper-distance`: upper triangular distance matrix; similar to the previous, but for the entries above the diagonal; suitable for output from the MATLAB functions `pdist` or `seqpdist`, exported to a CSV file. - `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 + - `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. - `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). + - `sparse`: sparse triplet format; a whitespace separated list of entries of a sparse distance matrix, one entry per line, each of the form *i j d(i,j)* specifying the distance between points *i* and *j*. Each pair of points should appear in the file at most once. - `--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. + - `--modulus p`: compute homology with coefficients in the prime field Z/*p*Z (only available when built with the option `USE_COEFFICIENTS`). - `--ratio `: only show persistence pairs with death/birth ratio > *r*. @@ -95,12 +96,10 @@ The following experimental features are currently available in separate branches - `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. - ### License -Ripser is licensed under the [MIT] license (`COPYING.txt`), with an extra clause (`CONTRIBUTING.txt`) clarifying the license for modifications released without an explicit written license agreement. Please contact the author if you want to use Ripser in your software under a different license. +Ripser is licensed under the [MIT] license (`COPYING.txt`), with an extra clause (`CONTRIBUTING.txt`) clarifying the license for modifications released without an explicit written license agreement. Please contact the author if you want to use Ripser in your software under a different license. [Ulrich Bauer]: [live.ripser.org]: diff --git a/ripser.cpp b/ripser.cpp index 8fb0cd1..c1d4fb3 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -998,17 +998,17 @@ void print_usage_and_exit(int exit_code) { << " distance (full distance matrix)" << std::endl << " point-cloud (point cloud in Euclidean space)" << std::endl << " dipha (distance matrix in DIPHA file format)" << std::endl - << " sparse (sparse distance matrix in Sparse Triplet format)" + << " sparse (sparse distance matrix in sparse triplet format)" << std::endl << " binary (lower triangular distance matrix in binary format)" << std::endl << " --dim compute persistent homology up to dimension k" << std::endl << " --threshold compute Rips complexes up to diameter t" << std::endl #ifdef USE_COEFFICIENTS - << " --modulus

compute homology with coefficients in the prime field Z/pZ" + << " --modulus

compute homology with coefficients in the prime field Z/pZ" << std::endl #endif - << 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 + << std::endl; exit(exit_code); } -- cgit v1.2.3 From 3fd98444e764f490752aac427011180c19e69e31 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Mon, 19 Aug 2019 12:42:19 +0200 Subject: fixed size_t / index_t conversion warnings --- ripser.cpp | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) (limited to 'ripser.cpp') diff --git a/ripser.cpp b/ripser.cpp index c1d4fb3..9418141 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -60,7 +60,6 @@ template class hash_map : public google::dense_hash_map { public: explicit hash_map() : google::dense_hash_map() { this->set_empty_key(-1); } - inline void reserve(size_t hint) { this->resize(hint); } }; #else @@ -221,8 +220,8 @@ template struct compressed_distance_matrix { : distances(mat.size() * (mat.size() - 1) / 2), rows(mat.size()) { init_rows(); - for (index_t i = 1; i < size(); ++i) - for (index_t j = 0; j < i; ++j) rows[i][j] = mat(i, j); + for (size_t i = 1; i < size(); ++i) + for (size_t j = 0; j < i; ++j) rows[i][j] = mat(i, j); } value_t operator()(const index_t i, const index_t j) const; @@ -235,7 +234,7 @@ typedef compressed_distance_matrix compressed_upper_distance_m template <> void compressed_lower_distance_matrix::init_rows() { value_t* pointer = &distances[0]; - for (index_t i = 1; i < size(); ++i) { + for (size_t i = 1; i < size(); ++i) { rows[i] = pointer; pointer += i; } @@ -243,7 +242,7 @@ template <> void compressed_lower_distance_matrix::init_rows() { template <> void compressed_upper_distance_matrix::init_rows() { value_t* pointer = &distances[0] - 1; - for (index_t i = 0; i < size() - 1; ++i) { + for (size_t i = 0; i < size() - 1; ++i) { rows[i] = pointer; pointer += size() - i - 2; } @@ -275,8 +274,8 @@ struct sparse_distance_matrix { sparse_distance_matrix(const DistanceMatrix& mat, const value_t threshold) : neighbors(mat.size()), num_edges(0) { - for (index_t i = 0; i < size(); ++i) - for (index_t j = 0; j < size(); ++j) + for (size_t i = 0; i < size(); ++i) + for (size_t j = 0; j < size(); ++j) if (i != j && mat(i, j) <= threshold) { ++num_edges; neighbors[i].push_back({j, mat(i, j)}); @@ -608,7 +607,7 @@ public: #ifdef INDICATE_PROGRESS std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step; #endif - for (index_t index_column_to_reduce = 0; index_column_to_reduce < columns_to_reduce.size(); + for (size_t index_column_to_reduce = 0; index_column_to_reduce < columns_to_reduce.size(); ++index_column_to_reduce) { diameter_entry_t column_to_reduce(columns_to_reduce[index_column_to_reduce], 1); @@ -889,7 +888,7 @@ sparse_distance_matrix read_sparse_distance_matrix(std::istream& input_stream) { } } - for (index_t i = 0; i < neighbors.size(); ++i) + for (size_t i = 0; i < neighbors.size(); ++i) std::sort(neighbors[i].begin(), neighbors[i].end()); return sparse_distance_matrix(std::move(neighbors), num_edges); @@ -1097,9 +1096,9 @@ int main(int argc, char** argv) { if (threshold == std::numeric_limits::max()) { value_t enclosing_radius = std::numeric_limits::infinity(); - for (index_t i = 0; i < dist.size(); ++i) { + for (size_t i = 0; i < dist.size(); ++i) { value_t r_i = -std::numeric_limits::infinity(); - for (index_t j = 0; j < dist.size(); ++j) r_i = std::max(r_i, dist(i, j)); + for (size_t j = 0; j < dist.size(); ++j) r_i = std::max(r_i, dist(i, j)); enclosing_radius = std::min(enclosing_radius, r_i); } threshold = enclosing_radius; -- cgit v1.2.3 From 3c6186b8450ef1c357b21042d9e963060a904edc Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Mon, 19 Aug 2019 12:49:33 +0200 Subject: code cleanup pop_pivot --- ripser.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) (limited to 'ripser.cpp') diff --git a/ripser.cpp b/ripser.cpp index 9418141..4d8e8f3 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -521,13 +521,13 @@ public: if (get_coefficient(pivot) == 0) pivot = column.top(); else if (get_index(column.top()) != get_index(pivot)) - break; + return pivot; else set_coefficient(pivot, (get_coefficient(pivot) + get_coefficient(column.top())) % modulus); column.pop(); } - if (get_coefficient(pivot) == 0) pivot = -1; + return (get_coefficient(pivot) == 0) ? -1 : pivot; #else while (!column.empty()) { pivot = column.top(); @@ -535,9 +535,8 @@ public: if (column.empty() || get_index(column.top()) != get_index(pivot)) return pivot; column.pop(); } - pivot = -1; + return -1; #endif - return pivot; } template diameter_entry_t get_pivot(Column& column) { -- cgit v1.2.3 From 94eb1f7db4421e324a827590f60aa236af56666c Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Tue, 20 Aug 2019 22:58:21 +0200 Subject: coboundary enumerator cleanup --- ripser.cpp | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) (limited to 'ripser.cpp') diff --git a/ripser.cpp b/ripser.cpp index 4d8e8f3..5549016 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -718,18 +718,17 @@ public: } bool has_next(bool all_cofacets = true) { - while ((v != -1) && (binomial_coeff(v, k) <= idx_below)) { - if (!all_cofacets) return false; + return (v >= k && (all_cofacets || binomial_coeff(v, k) > idx_below)); + } + + diameter_entry_t next() { + while ((binomial_coeff(v, k) <= idx_below)) { idx_below -= binomial_coeff(v, k); idx_above += binomial_coeff(v, k + 1); --v; --k; assert(k != -1); } - return v != -1; - } - - diameter_entry_t next() { value_t cofacet_diameter = get_diameter(simplex); for (index_t w : vertices) cofacet_diameter = std::max(cofacet_diameter, dist(v, w)); index_t cofacet_index = idx_above + binomial_coeff(v--, k + 1) + idx_below; @@ -754,10 +753,10 @@ template <> class ripser::simplex_coboundary_enumerator public: 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), - k(_dim + 1), - vertices(_dim + 1), simplex(_simplex), modulus(parent.modulus), dist(parent.dist), binomial_coeff(parent.binomial_coeff), - neighbor_it(dist.neighbor_it), neighbor_end(dist.neighbor_end) { + : parent(_parent), idx_below(get_index(_simplex)), idx_above(0), k(_dim + 1), + vertices(_dim + 1), simplex(_simplex), modulus(parent.modulus), dist(parent.dist), + binomial_coeff(parent.binomial_coeff), neighbor_it(dist.neighbor_it), + neighbor_end(dist.neighbor_end) { neighbor_it.clear(); neighbor_end.clear(); @@ -1003,7 +1002,8 @@ void print_usage_and_exit(int exit_code) { << " --dim compute persistent homology up to dimension k" << std::endl << " --threshold compute Rips complexes up to diameter t" << std::endl #ifdef USE_COEFFICIENTS - << " --modulus

compute homology with coefficients in the prime field Z/pZ" << std::endl + << " --modulus

compute homology with coefficients in the prime field Z/pZ" + << std::endl #endif << " --ratio only show persistence pairs with death/birth ratio > r" << std::endl << std::endl; -- cgit v1.2.3