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