From 11dfa016218793929b6a8b7f8c93a3669f3229d5 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 29 Dec 2016 17:52:49 +0100 Subject: code reformatting --- .clang-format | 2 +- ripser.cpp | 469 ++++++++++++++++++++++++++++++++-------------------------- 2 files changed, 259 insertions(+), 212 deletions(-) diff --git a/.clang-format b/.clang-format index 598d7c2..1975b19 100644 --- a/.clang-format +++ b/.clang-format @@ -1,7 +1,7 @@ BasedOnStyle: LLVM IndentWidth: 4 TabWidth: 4 -ColumnLimit: 120 +ColumnLimit: 100 AccessModifierOffset: -4 AllowShortIfStatementsOnASingleLine: true AllowShortLoopsOnASingleLine: true diff --git a/ripser.cpp b/ripser.cpp index 052c07d..02bfac9 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -52,6 +52,7 @@ typedef int16_t coefficient_t; class binomial_coeff_table { std::vector> B; + public: binomial_coeff_table(index_t n, index_t k) : B(n + 1) { for (index_t i = 0; i <= n; i++) { @@ -87,19 +88,21 @@ std::vector multiplicative_inverse_vector(const coefficient_t m) return inverse; } - #ifdef USE_COEFFICIENTS struct __attribute__((packed)) entry_t { index_t index : 8 * (sizeof(index_t) - sizeof(coefficient_t)); coefficient_t coefficient; - entry_t(index_t _index, coefficient_t _coefficient) : index(_index), coefficient(_coefficient) {} + entry_t(index_t _index, coefficient_t _coefficient) + : index(_index), coefficient(_coefficient) {} entry_t(index_t _index) : index(_index), coefficient(1) {} entry_t() : index(0), coefficient(1) {} }; 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); } +entry_t make_entry(index_t _index, coefficient_t _coefficient) { + return entry_t(_index, _coefficient); +} index_t get_index(entry_t e) { return e.index; } index_t get_coefficient(entry_t e) { return e.coefficient; } void set_coefficient(entry_t& e, const coefficient_t c) { e.coefficient = c; } @@ -149,9 +152,13 @@ public: const entry_t& get_entry(const diameter_entry_t& p) { return p.second; } entry_t& get_entry(diameter_entry_t& p) { return p.second; } const index_t get_index(const diameter_entry_t& p) { return get_index(get_entry(p)); } -const coefficient_t get_coefficient(const diameter_entry_t& p) { return get_coefficient(get_entry(p)); } +const coefficient_t get_coefficient(const diameter_entry_t& p) { + return get_coefficient(get_entry(p)); +} const value_t& get_diameter(const diameter_entry_t& p) { return p.first; } -void set_coefficient(diameter_entry_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); } +void set_coefficient(diameter_entry_t& p, const coefficient_t c) { + set_coefficient(get_entry(p), c); +} template struct greater_diameter_or_smaller_index { bool operator()(const Entry& a, const Entry& b) { @@ -205,11 +212,13 @@ template <> void compressed_distance_matrix::init_rows() { } } -template <> value_t compressed_distance_matrix::operator()(index_t i, index_t j) const { +template <> +value_t compressed_distance_matrix::operator()(index_t i, index_t j) const { return i == j ? 0 : i > j ? rows[j][i] : rows[i][j]; } -template <> value_t compressed_distance_matrix::operator()(index_t i, index_t j) const { +template <> +value_t compressed_distance_matrix::operator()(index_t i, index_t j) const { return i == j ? 0 : i < j ? rows[j][i] : rows[i][j]; } @@ -220,12 +229,13 @@ class euclidean_distance_matrix { public: std::vector> points; - euclidean_distance_matrix(std::vector>&& _points) : points(std::move(_points)) {} + euclidean_distance_matrix(std::vector>&& _points) + : points(std::move(_points)) {} value_t operator()(const index_t i, const index_t j) const { - return std::sqrt(std::inner_product(points[i].begin(), points[i].end(), points[j].begin(), value_t(), - std::plus(), - [](value_t u, value_t v) { return (u - v) * (u - v); })); + return std::sqrt(std::inner_product( + points[i].begin(), points[i].end(), points[j].begin(), value_t(), std::plus(), + [](value_t u, value_t v) { return (u - v) * (u - v); })); } size_t size() const { return points.size(); } @@ -350,7 +360,8 @@ public: } }; -template void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) { +template +void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) { entry_t e = make_entry(i, c); column.push(std::make_pair(diameter, e)); } @@ -365,11 +376,12 @@ class ripser { mutable std::vector vertices; public: - ripser(compressed_lower_distance_matrix&& _dist, index_t _dim_max, value_t _threshold, coefficient_t _modulus) - : dist(std::move(_dist)), dim_max(std::min(_dim_max, index_t(dist.size() - 2))), n(dist.size()), threshold(_threshold), - modulus(_modulus), binomial_coeff(n, dim_max + 2), + ripser(compressed_lower_distance_matrix&& _dist, index_t _dim_max, value_t _threshold, + coefficient_t _modulus) + : dist(std::move(_dist)), dim_max(std::min(_dim_max, index_t(dist.size() - 2))), + n(dist.size()), threshold(_threshold), 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 { if (binomial_coeff(v, k) > idx) { index_t count = v; @@ -387,10 +399,10 @@ public: assert(binomial_coeff(v, k) <= idx && binomial_coeff(v + 1, k) > idx); return v; } - + template OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v, - OutputIterator out) const { + OutputIterator out) const { --v; for (index_t k = dim + 1; k > 0; --k) { get_next_vertex(v, idx, k); @@ -399,7 +411,7 @@ public: } return out; } - + value_t compute_diameter(const index_t index, index_t dim) const { value_t diam = -std::numeric_limits::infinity(); @@ -407,10 +419,12 @@ public: get_simplex_vertices(index, dim, dist.size(), std::back_inserter(vertices)); for (index_t i = 0; i <= dim; ++i) - for (index_t j = 0; j < i; ++j) { diam = std::max(diam, dist(vertices[i], vertices[j])); } + for (index_t j = 0; j < i; ++j) { + diam = std::max(diam, dist(vertices[i], vertices[j])); + } return diam; } - + class simplex_coboundary_enumerator { private: index_t idx_below, idx_above, v, k; @@ -419,14 +433,16 @@ public: const coefficient_t modulus; const compressed_lower_distance_matrix& dist; const binomial_coeff_table& binomial_coeff; - + public: - simplex_coboundary_enumerator(const diameter_entry_t _simplex, 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), - binomial_coeff(parent.binomial_coeff) { + simplex_coboundary_enumerator(const diameter_entry_t _simplex, 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), + binomial_coeff(parent.binomial_coeff) { parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.begin()); } - + bool has_next() { while ((v != -1) && (binomial_coeff(v, k) <= idx_below)) { idx_below -= binomial_coeff(v, k); @@ -437,276 +453,248 @@ public: } return v != -1; } - + 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 = (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; + coefficient_t coface_coefficient = + (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; return diameter_entry_t(coface_diameter, coface_index, coface_coefficient); } }; - void assemble_columns_to_reduce(std::vector& columns_to_reduce, - hash_map& pivot_column_index, index_t dim) { - index_t num_simplices = binomial_coeff(n, dim + 1); + hash_map& pivot_column_index, index_t dim); + + void compute_pairs(std::vector& columns_to_reduce, + hash_map& pivot_column_index, index_t dim); + + void compute_barcodes(); +}; - columns_to_reduce.clear(); +void ripser::assemble_columns_to_reduce(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"; + 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)); + 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) % 1000 == 0) - std::cout << "\033[K" - << "assembled " << columns_to_reduce.size() << " out of " << (index + 1) << "/" - << num_simplices << " columns" << std::flush << "\r"; + if ((index + 1) % 1000 == 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"; + 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()); + std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), + greater_diameter_or_smaller_index()); #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif - } +} - void compute_pairs(std::vector& columns_to_reduce, hash_map& pivot_column_index, - index_t dim) { +void ripser::compute_pairs(std::vector& columns_to_reduce, + hash_map& pivot_column_index, index_t dim) { #ifdef PRINT_PERSISTENCE_PAIRS - std::cout << "persistence intervals in dim " << dim << ":" << std::endl; + std::cout << "persistence intervals in dim " << dim << ":" << std::endl; #endif #ifdef ASSEMBLE_REDUCTION_MATRIX - compressed_sparse_matrix reduction_coefficients; + compressed_sparse_matrix reduction_coefficients; #else #ifdef USE_COEFFICIENTS - std::vector reduction_coefficients; + std::vector reduction_coefficients; #endif #endif - std::vector coface_entries; + std::vector coface_entries; - for (index_t i = 0; i < columns_to_reduce.size(); ++i) { - auto column_to_reduce = columns_to_reduce[i]; + for (index_t i = 0; i < columns_to_reduce.size(); ++i) { + auto column_to_reduce = columns_to_reduce[i]; #ifdef ASSEMBLE_REDUCTION_MATRIX - std::priority_queue, - greater_diameter_or_smaller_index> - reduction_column; + std::priority_queue, + greater_diameter_or_smaller_index> + reduction_column; #endif - std::priority_queue, - greater_diameter_or_smaller_index> - working_coboundary; + std::priority_queue, + greater_diameter_or_smaller_index> + working_coboundary; - value_t diameter = get_diameter(column_to_reduce); + value_t diameter = get_diameter(column_to_reduce); #ifdef INDICATE_PROGRESS - if ((i + 1) % 1000 == 0) - std::cout << "\033[K" - << "reducing column " << i + 1 << "/" << columns_to_reduce.size() << " (diameter " << diameter - << ")" << std::flush << "\r"; + if ((i + 1) % 1000 == 0) + std::cout << "\033[K" + << "reducing column " << i + 1 << "/" << columns_to_reduce.size() + << " (diameter " << diameter << ")" << std::flush << "\r"; #endif - index_t j = i; + index_t j = i; - // start with a dummy pivot entry with coefficient -1 in order to initialize - // working_coboundary with the coboundary of the simplex with index column_to_reduce - diameter_entry_t pivot(0, -1, -1 + modulus); + // start with a dummy pivot entry with coefficient -1 in order to initialize + // working_coboundary with the coboundary of the simplex with index column_to_reduce + diameter_entry_t pivot(0, -1, -1 + modulus); #ifdef ASSEMBLE_REDUCTION_MATRIX - // initialize reduction_coefficients as identity matrix - reduction_coefficients.append_column(); + // initialize reduction_coefficients as identity matrix + reduction_coefficients.append_column(); #endif #ifdef USE_COEFFICIENTS - reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, 1)); + reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, 1)); #endif - bool might_be_apparent_pair = (i == j); + bool might_be_apparent_pair = (i == j); - do { - const coefficient_t factor = modulus - get_coefficient(pivot); + do { + const coefficient_t factor = modulus - get_coefficient(pivot); #ifdef ASSEMBLE_REDUCTION_MATRIX #ifdef USE_COEFFICIENTS - auto coeffs_begin = reduction_coefficients.cbegin(j), coeffs_end = reduction_coefficients.cend(j); + auto coeffs_begin = reduction_coefficients.cbegin(j), + coeffs_end = reduction_coefficients.cend(j); #else - std::vector coeffs; - coeffs.push_back(columns_to_reduce[j]); - for (auto it = reduction_coefficients.cbegin(j); it != reduction_coefficients.cend(j); ++it) - coeffs.push_back(*it); - auto coeffs_begin = coeffs.begin(), coeffs_end = coeffs.end(); + std::vector coeffs; + coeffs.push_back(columns_to_reduce[j]); + for (auto it = reduction_coefficients.cbegin(j); it != reduction_coefficients.cend(j); + ++it) + coeffs.push_back(*it); + auto coeffs_begin = coeffs.begin(), coeffs_end = coeffs.end(); #endif #else #ifdef USE_COEFFICIENTS - auto coeffs_begin = &reduction_coefficients[j], coeffs_end = &reduction_coefficients[j] + 1; + auto coeffs_begin = &reduction_coefficients[j], + coeffs_end = &reduction_coefficients[j] + 1; #else - auto coeffs_begin = &columns_to_reduce[j], coeffs_end = &columns_to_reduce[j] + 1; + auto coeffs_begin = &columns_to_reduce[j], coeffs_end = &columns_to_reduce[j] + 1; #endif #endif - for (auto it = coeffs_begin; it != coeffs_end; ++it) { - diameter_entry_t simplex = *it; - set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); + for (auto it = coeffs_begin; it != coeffs_end; ++it) { + diameter_entry_t simplex = *it; + set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); #ifdef ASSEMBLE_REDUCTION_MATRIX - reduction_column.push(simplex); -#endif - - 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()) { - pivot = coface; - goto found_persistence_pair; - } - might_be_apparent_pair = false; + reduction_column.push(simplex); +#endif + + 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()) { + pivot = coface; + goto found_persistence_pair; } + might_be_apparent_pair = false; } } - for (auto coface : coface_entries) working_coboundary.push(coface); } + for (auto coface : coface_entries) working_coboundary.push(coface); + } - pivot = get_pivot(working_coboundary, modulus); + pivot = get_pivot(working_coboundary, modulus); - if (get_index(pivot) != -1) { - auto pair = pivot_column_index.find(get_index(pivot)); + if (get_index(pivot) != -1) { + auto pair = pivot_column_index.find(get_index(pivot)); - if (pair != pivot_column_index.end()) { - j = pair->second; - continue; - } - } else { + if (pair != pivot_column_index.end()) { + j = pair->second; + continue; + } + } else { #ifdef PRINT_PERSISTENCE_PAIRS #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif - std::cout << " [" << diameter << ", )" << std::endl << std::flush; + std::cout << " [" << diameter << ", )" << std::endl << std::flush; #endif - break; - } + break; + } - found_persistence_pair: + found_persistence_pair: #ifdef PRINT_PERSISTENCE_PAIRS - value_t death = get_diameter(pivot); - if (diameter != death) { + value_t death = get_diameter(pivot); + if (diameter != death) { #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif - std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush; - } + std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush; + } #endif - pivot_column_index.insert(std::make_pair(get_index(pivot), i)); + pivot_column_index.insert(std::make_pair(get_index(pivot), i)); #ifdef USE_COEFFICIENTS - const coefficient_t inverse = multiplicative_inverse[get_coefficient(pivot)]; + const coefficient_t inverse = multiplicative_inverse[get_coefficient(pivot)]; #endif #ifdef ASSEMBLE_REDUCTION_MATRIX // replace current column of reduction_coefficients (with a single diagonal 1 entry) // by reduction_column (possibly with a different entry on the diagonal) #ifdef USE_COEFFICIENTS - reduction_coefficients.pop_back(); + reduction_coefficients.pop_back(); #else - pop_pivot(reduction_column, modulus); + pop_pivot(reduction_column, modulus); #endif - while (true) { - diameter_entry_t e = pop_pivot(reduction_column, modulus); - if (get_index(e) == -1) break; + while (true) { + diameter_entry_t e = pop_pivot(reduction_column, modulus); + if (get_index(e) == -1) break; #ifdef USE_COEFFICIENTS - set_coefficient(e, inverse * get_coefficient(e) % modulus); - assert(get_coefficient(e) > 0); + set_coefficient(e, inverse * get_coefficient(e) % modulus); + assert(get_coefficient(e) > 0); #endif - reduction_coefficients.push_back(e); - } + reduction_coefficients.push_back(e); + } #else #ifdef USE_COEFFICIENTS - reduction_coefficients.pop_back(); - reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, inverse)); + reduction_coefficients.pop_back(); + reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, inverse)); #endif #endif - break; - } while (true); - } - -#ifdef INDICATE_PROGRESS - std::cout << "\033[K"; -#endif + break; + } while (true); } - void compute_barcodes() { - - std::vector columns_to_reduce; - - { - union_find dset(n); - std::vector edges; - 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)); - } - 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) { - vertices_of_edge.clear(); - get_simplex_vertices(get_index(e), 1, n, std::back_inserter(vertices_of_edge)); - index_t u = dset.find(vertices_of_edge[0]), v = dset.find(vertices_of_edge[1]); - - if (u != v) { -#ifdef PRINT_PERSISTENCE_PAIRS - if (get_diameter(e) != 0) std::cout << " [0," << get_diameter(e) << ")" << std::endl; -#endif - dset.link(u, v); - } else - columns_to_reduce.push_back(e); - } - std::reverse(columns_to_reduce.begin(), columns_to_reduce.end()); - -#ifdef PRINT_PERSISTENCE_PAIRS - for (index_t i = 0; i < n; ++i) - if (dset.find(i) == i) std::cout << " [0, )" << std::endl << std::flush; +#ifdef INDICATE_PROGRESS + std::cout << "\033[K"; #endif - } - - for (index_t dim = 1; dim <= dim_max; ++dim) { - hash_map pivot_column_index; - pivot_column_index.reserve(columns_to_reduce.size()); - - compute_pairs(columns_to_reduce, pivot_column_index, dim); +} - if (dim < dim_max) { assemble_columns_to_reduce(columns_to_reduce, pivot_column_index, dim + 1); } - } - } +enum file_format { + LOWER_DISTANCE_MATRIX, + UPPER_DISTANCE_MATRIX, + DISTANCE_MATRIX, + POINT_CLOUD, + DIPHA, + RIPSER }; -enum file_format { LOWER_DISTANCE_MATRIX, UPPER_DISTANCE_MATRIX, DISTANCE_MATRIX, POINT_CLOUD, DIPHA, RIPSER }; - template T read(std::istream& s) { T result; s.read(reinterpret_cast(&result), sizeof(T)); @@ -733,7 +721,8 @@ compressed_lower_distance_matrix read_point_cloud(std::istream& input_stream) { index_t n = eucl_dist.size(); - std::cout << "point cloud with " << n << " points in dimension " << eucl_dist.points.front().size() << std::endl; + std::cout << "point cloud with " << n << " points in dimension " + << eucl_dist.points.front().size() << std::endl; std::vector distances; @@ -830,26 +819,30 @@ compressed_lower_distance_matrix read_file(std::istream& input_stream, file_form } void print_usage_and_exit(int exit_code) { - std::cerr << "Usage: " - << "ripser " - << "[options] [filename]" << std::endl - << std::endl - << "Options:" << std::endl - << std::endl - << " --help print this screen" << std::endl - << " --format use the specified file format for the input. Options are:" << std::endl - << " lower-distance (lower triangular distance matrix; default)" << std::endl - << " upper-distance (upper triangular distance matrix)" << std::endl - << " distance (full distance matrix)" << std::endl - << " point-cloud (point cloud in Euclidean space)" << std::endl - << " dipha (distance matrix in DIPHA file format)" << std::endl - << " ripser (distance matrix in Ripser binary file format)" << std::endl - << " --dim compute persistent homology up to dimension " << std::endl - << " --threshold compute Rips complexes up to diameter " << std::endl + std::cerr + << "Usage: " + << "ripser " + << "[options] [filename]" << std::endl + << std::endl + << "Options:" << std::endl + << std::endl + << " --help print this screen" << std::endl + << " --format use the specified file format for the input. Options are:" + << std::endl + << " lower-distance (lower triangular distance matrix; default)" + << std::endl + << " upper-distance (upper triangular distance matrix)" << std::endl + << " distance (full distance matrix)" << std::endl + << " point-cloud (point cloud in Euclidean space)" << std::endl + << " dipha (distance matrix in DIPHA file format)" << std::endl + << " ripser (distance matrix in Ripser binary file format)" + << std::endl + << " --dim compute persistent homology up to dimension " << std::endl + << " --threshold compute Rips complexes up to diameter " << std::endl #ifdef USE_COEFFICIENTS - << " --modulus

compute homology with coefficients in the prime field Z/

Z" + << " --modulus

compute homology with coefficients in the prime field Z/

Z" #endif - << std::endl; + << std::endl; exit(exit_code); } @@ -923,7 +916,61 @@ int main(int argc, char** argv) { std::cout << "distance matrix with " << dist.size() << " points" << std::endl; auto value_range = std::minmax_element(dist.distances.begin(), dist.distances.end()); - std::cout << "value range: [" << *value_range.first << "," << *value_range.second << "]" << std::endl; + std::cout << "value range: [" << *value_range.first << "," << *value_range.second << "]" + << std::endl; ripser(std::move(dist), dim_max, threshold, modulus).compute_barcodes(); } + +void ripser::compute_barcodes() { + + std::vector columns_to_reduce; + + { + union_find dset(n); + std::vector edges; + 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)); + } + 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) { + vertices_of_edge.clear(); + get_simplex_vertices(get_index(e), 1, n, std::back_inserter(vertices_of_edge)); + index_t u = dset.find(vertices_of_edge[0]), v = dset.find(vertices_of_edge[1]); + + if (u != v) { +#ifdef PRINT_PERSISTENCE_PAIRS + if (get_diameter(e) != 0) + std::cout << " [0," << get_diameter(e) << ")" << std::endl; +#endif + dset.link(u, v); + } else + columns_to_reduce.push_back(e); + } + std::reverse(columns_to_reduce.begin(), columns_to_reduce.end()); + +#ifdef PRINT_PERSISTENCE_PAIRS + for (index_t i = 0; i < n; ++i) + if (dset.find(i) == i) std::cout << " [0, )" << std::endl << std::flush; +#endif + } + + for (index_t dim = 1; dim <= dim_max; ++dim) { + hash_map pivot_column_index; + pivot_column_index.reserve(columns_to_reduce.size()); + + compute_pairs(columns_to_reduce, pivot_column_index, dim); + + if (dim < dim_max) { + assemble_columns_to_reduce(columns_to_reduce, pivot_column_index, dim + 1); + } + } +} -- cgit v1.2.3 From 2bed1bcf7f3d654b0f8332e003bb3ec602b84681 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 29 Dec 2016 17:59:03 +0100 Subject: code reformatting --- ripser.cpp | 247 ++++++++++++++++++++++++++++++------------------------------- 1 file changed, 121 insertions(+), 126 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 02bfac9..e17bc6a 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -464,227 +464,222 @@ public: } }; - void assemble_columns_to_reduce(std::vector& columns_to_reduce, - hash_map& pivot_column_index, index_t dim); - - void compute_pairs(std::vector& columns_to_reduce, - hash_map& pivot_column_index, index_t dim); - void compute_barcodes(); -}; -void ripser::assemble_columns_to_reduce(std::vector& columns_to_reduce, - hash_map& pivot_column_index, - index_t dim) { - index_t num_simplices = binomial_coeff(n, dim + 1); + void assemble_columns_to_reduce(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(); + columns_to_reduce.clear(); #ifdef INDICATE_PROGRESS - std::cout << "\033[K" - << "assembling " << num_simplices << " columns" << std::flush << "\r"; + 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)); + 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) % 1000 == 0) - std::cout << "\033[K" - << "assembled " << columns_to_reduce.size() << " out of " << (index + 1) - << "/" << num_simplices << " columns" << std::flush << "\r"; + if ((index + 1) % 1000 == 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"; + 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()); + std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), + greater_diameter_or_smaller_index()); #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif -} + } -void ripser::compute_pairs(std::vector& columns_to_reduce, - hash_map& pivot_column_index, index_t dim) { + void compute_pairs(std::vector& columns_to_reduce, + hash_map& pivot_column_index, index_t dim) { #ifdef PRINT_PERSISTENCE_PAIRS - std::cout << "persistence intervals in dim " << dim << ":" << std::endl; + std::cout << "persistence intervals in dim " << dim << ":" << std::endl; #endif #ifdef ASSEMBLE_REDUCTION_MATRIX - compressed_sparse_matrix reduction_coefficients; + compressed_sparse_matrix reduction_coefficients; #else #ifdef USE_COEFFICIENTS - std::vector reduction_coefficients; + std::vector reduction_coefficients; #endif #endif - std::vector coface_entries; + std::vector coface_entries; - for (index_t i = 0; i < columns_to_reduce.size(); ++i) { - auto column_to_reduce = columns_to_reduce[i]; + for (index_t i = 0; i < columns_to_reduce.size(); ++i) { + auto column_to_reduce = columns_to_reduce[i]; #ifdef ASSEMBLE_REDUCTION_MATRIX - std::priority_queue, - greater_diameter_or_smaller_index> - reduction_column; + std::priority_queue, + greater_diameter_or_smaller_index> + reduction_column; #endif - std::priority_queue, - greater_diameter_or_smaller_index> - working_coboundary; + std::priority_queue, + greater_diameter_or_smaller_index> + working_coboundary; - value_t diameter = get_diameter(column_to_reduce); + value_t diameter = get_diameter(column_to_reduce); #ifdef INDICATE_PROGRESS - if ((i + 1) % 1000 == 0) - std::cout << "\033[K" - << "reducing column " << i + 1 << "/" << columns_to_reduce.size() - << " (diameter " << diameter << ")" << std::flush << "\r"; + if ((i + 1) % 1000 == 0) + std::cout << "\033[K" + << "reducing column " << i + 1 << "/" << columns_to_reduce.size() + << " (diameter " << diameter << ")" << std::flush << "\r"; #endif - index_t j = i; + index_t j = i; - // start with a dummy pivot entry with coefficient -1 in order to initialize - // working_coboundary with the coboundary of the simplex with index column_to_reduce - diameter_entry_t pivot(0, -1, -1 + modulus); + // start with a dummy pivot entry with coefficient -1 in order to initialize + // working_coboundary with the coboundary of the simplex with index column_to_reduce + diameter_entry_t pivot(0, -1, -1 + modulus); #ifdef ASSEMBLE_REDUCTION_MATRIX - // initialize reduction_coefficients as identity matrix - reduction_coefficients.append_column(); + // initialize reduction_coefficients as identity matrix + reduction_coefficients.append_column(); #endif #ifdef USE_COEFFICIENTS - reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, 1)); + reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, 1)); #endif - bool might_be_apparent_pair = (i == j); + bool might_be_apparent_pair = (i == j); - do { - const coefficient_t factor = modulus - get_coefficient(pivot); + do { + const coefficient_t factor = modulus - get_coefficient(pivot); #ifdef ASSEMBLE_REDUCTION_MATRIX #ifdef USE_COEFFICIENTS - auto coeffs_begin = reduction_coefficients.cbegin(j), - coeffs_end = reduction_coefficients.cend(j); + auto coeffs_begin = reduction_coefficients.cbegin(j), + coeffs_end = reduction_coefficients.cend(j); #else - std::vector coeffs; - coeffs.push_back(columns_to_reduce[j]); - for (auto it = reduction_coefficients.cbegin(j); it != reduction_coefficients.cend(j); - ++it) - coeffs.push_back(*it); - auto coeffs_begin = coeffs.begin(), coeffs_end = coeffs.end(); + std::vector coeffs; + coeffs.push_back(columns_to_reduce[j]); + for (auto it = reduction_coefficients.cbegin(j); + it != reduction_coefficients.cend(j); ++it) + coeffs.push_back(*it); + auto coeffs_begin = coeffs.begin(), coeffs_end = coeffs.end(); #endif #else #ifdef USE_COEFFICIENTS - auto coeffs_begin = &reduction_coefficients[j], - coeffs_end = &reduction_coefficients[j] + 1; + auto coeffs_begin = &reduction_coefficients[j], + coeffs_end = &reduction_coefficients[j] + 1; #else - auto coeffs_begin = &columns_to_reduce[j], coeffs_end = &columns_to_reduce[j] + 1; + auto coeffs_begin = &columns_to_reduce[j], coeffs_end = &columns_to_reduce[j] + 1; #endif #endif - for (auto it = coeffs_begin; it != coeffs_end; ++it) { - diameter_entry_t simplex = *it; - set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); + for (auto it = coeffs_begin; it != coeffs_end; ++it) { + diameter_entry_t simplex = *it; + set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); #ifdef ASSEMBLE_REDUCTION_MATRIX - reduction_column.push(simplex); -#endif - - 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()) { - pivot = coface; - goto found_persistence_pair; + reduction_column.push(simplex); +#endif + + 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()) { + pivot = coface; + goto found_persistence_pair; + } + might_be_apparent_pair = false; } - might_be_apparent_pair = false; } } + for (auto coface : coface_entries) working_coboundary.push(coface); } - for (auto coface : coface_entries) working_coboundary.push(coface); - } - pivot = get_pivot(working_coboundary, modulus); + pivot = get_pivot(working_coboundary, modulus); - if (get_index(pivot) != -1) { - auto pair = pivot_column_index.find(get_index(pivot)); + if (get_index(pivot) != -1) { + auto pair = pivot_column_index.find(get_index(pivot)); - if (pair != pivot_column_index.end()) { - j = pair->second; - continue; - } - } else { + if (pair != pivot_column_index.end()) { + j = pair->second; + continue; + } + } else { #ifdef PRINT_PERSISTENCE_PAIRS #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif - std::cout << " [" << diameter << ", )" << std::endl << std::flush; + std::cout << " [" << diameter << ", )" << std::endl << std::flush; #endif - break; - } + break; + } - found_persistence_pair: + found_persistence_pair: #ifdef PRINT_PERSISTENCE_PAIRS - value_t death = get_diameter(pivot); - if (diameter != death) { + value_t death = get_diameter(pivot); + if (diameter != death) { #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif - std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush; - } + std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush; + } #endif - pivot_column_index.insert(std::make_pair(get_index(pivot), i)); + pivot_column_index.insert(std::make_pair(get_index(pivot), i)); #ifdef USE_COEFFICIENTS - const coefficient_t inverse = multiplicative_inverse[get_coefficient(pivot)]; + const coefficient_t inverse = multiplicative_inverse[get_coefficient(pivot)]; #endif #ifdef ASSEMBLE_REDUCTION_MATRIX // replace current column of reduction_coefficients (with a single diagonal 1 entry) // by reduction_column (possibly with a different entry on the diagonal) #ifdef USE_COEFFICIENTS - reduction_coefficients.pop_back(); + reduction_coefficients.pop_back(); #else - pop_pivot(reduction_column, modulus); + pop_pivot(reduction_column, modulus); #endif - while (true) { - diameter_entry_t e = pop_pivot(reduction_column, modulus); - if (get_index(e) == -1) break; + while (true) { + diameter_entry_t e = pop_pivot(reduction_column, modulus); + if (get_index(e) == -1) break; #ifdef USE_COEFFICIENTS - set_coefficient(e, inverse * get_coefficient(e) % modulus); - assert(get_coefficient(e) > 0); + set_coefficient(e, inverse * get_coefficient(e) % modulus); + assert(get_coefficient(e) > 0); #endif - reduction_coefficients.push_back(e); - } + reduction_coefficients.push_back(e); + } #else #ifdef USE_COEFFICIENTS - reduction_coefficients.pop_back(); - reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, inverse)); + reduction_coefficients.pop_back(); + reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, inverse)); #endif #endif - break; - } while (true); - } + break; + } while (true); + } #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif -} + } +}; enum file_format { LOWER_DISTANCE_MATRIX, -- cgit v1.2.3