From 7d662086064455d2b409d3b5651a7ce009b0a6d4 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 29 Apr 2018 21:49:25 +0200 Subject: xcode project update --- ripser.xcodeproj/project.pbxproj | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/ripser.xcodeproj/project.pbxproj b/ripser.xcodeproj/project.pbxproj index 24c7d8b..b5f6881 100644 --- a/ripser.xcodeproj/project.pbxproj +++ b/ripser.xcodeproj/project.pbxproj @@ -88,7 +88,7 @@ 551018401BD63CB300990BFF /* Project object */ = { isa = PBXProject; attributes = { - LastUpgradeCheck = 0800; + LastUpgradeCheck = 0910; ORGANIZATIONNAME = "ulrich-bauer.org"; TargetAttributes = { 551018471BD63CB300990BFF = { @@ -133,14 +133,20 @@ CLANG_CXX_LIBRARY = "libc++"; CLANG_ENABLE_MODULES = YES; CLANG_ENABLE_OBJC_ARC = YES; + CLANG_WARN_BLOCK_CAPTURE_AUTORELEASING = YES; CLANG_WARN_BOOL_CONVERSION = YES; + CLANG_WARN_COMMA = YES; CLANG_WARN_CONSTANT_CONVERSION = YES; CLANG_WARN_DIRECT_OBJC_ISA_USAGE = YES_ERROR; CLANG_WARN_EMPTY_BODY = YES; CLANG_WARN_ENUM_CONVERSION = YES; CLANG_WARN_INFINITE_RECURSION = YES; CLANG_WARN_INT_CONVERSION = YES; + CLANG_WARN_NON_LITERAL_NULL_CONVERSION = YES; + CLANG_WARN_OBJC_LITERAL_CONVERSION = YES; CLANG_WARN_OBJC_ROOT_CLASS = YES_ERROR; + CLANG_WARN_RANGE_LOOP_ANALYSIS = YES; + CLANG_WARN_STRICT_PROTOTYPES = YES; CLANG_WARN_SUSPICIOUS_MOVE = YES; CLANG_WARN_UNREACHABLE_CODE = YES; CLANG_WARN__DUPLICATE_METHOD_MATCH = YES; @@ -178,14 +184,20 @@ CLANG_CXX_LIBRARY = "libc++"; CLANG_ENABLE_MODULES = YES; CLANG_ENABLE_OBJC_ARC = YES; + CLANG_WARN_BLOCK_CAPTURE_AUTORELEASING = YES; CLANG_WARN_BOOL_CONVERSION = YES; + CLANG_WARN_COMMA = YES; CLANG_WARN_CONSTANT_CONVERSION = YES; CLANG_WARN_DIRECT_OBJC_ISA_USAGE = YES_ERROR; CLANG_WARN_EMPTY_BODY = YES; CLANG_WARN_ENUM_CONVERSION = YES; CLANG_WARN_INFINITE_RECURSION = YES; CLANG_WARN_INT_CONVERSION = YES; + CLANG_WARN_NON_LITERAL_NULL_CONVERSION = YES; + CLANG_WARN_OBJC_LITERAL_CONVERSION = YES; CLANG_WARN_OBJC_ROOT_CLASS = YES_ERROR; + CLANG_WARN_RANGE_LOOP_ANALYSIS = YES; + CLANG_WARN_STRICT_PROTOTYPES = YES; CLANG_WARN_SUSPICIOUS_MOVE = YES; CLANG_WARN_UNREACHABLE_CODE = YES; CLANG_WARN__DUPLICATE_METHOD_MATCH = YES; -- cgit v1.2.3 From e10038a0fbd1bf2f8c9eec129b2d523d95a4ad0f Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 16 May 2018 09:25:46 +0200 Subject: simplify union-find --- ripser.cpp | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 6bd725e..bd1a7db 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -252,23 +252,16 @@ public: } index_t find(index_t x) { - index_t y = x, z = parent[y]; - while (z != y) { - y = z; - z = parent[y]; - } - y = parent[x]; - while (z != y) { - parent[x] = z; - x = y; - y = parent[x]; + index_t y = x, z; + while ((z = parent[y]) != y) y = z; + while ((z = parent[x]) != y) { + parent[x] = y; + x = z; } return z; } void link(index_t x, index_t y) { - x = find(x); - y = find(y); - if (x == y) return; + if ((x = find(x)) == (y = find(y))) return; if (rank[x] > rank[y]) parent[y] = x; else { -- cgit v1.2.3 From 932cc3028fcde38e6e4d6588f1f7bde76a344eb3 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sat, 19 May 2018 12:58:37 -0400 Subject: consolidation of branches --- ripser.cpp | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index bd1a7db..8f2d509 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -208,12 +208,14 @@ template <> void compressed_distance_matrix::init_rows() { } template <> -value_t compressed_distance_matrix::operator()(index_t i, index_t j) const { +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()(index_t i, index_t j) const { +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]; } @@ -226,10 +228,8 @@ public: euclidean_distance_matrix(std::vector>&& _points) : points(std::move(_points)) { - for (auto p: points) { - assert(p.size() == points.front().size()); - } - } + for (auto p : points) { assert(p.size() == points.front().size()); } + } value_t operator()(const index_t i, const index_t j) const { assert(i < points.size()); @@ -356,18 +356,20 @@ public: class ripser { compressed_lower_distance_matrix dist; - index_t dim_max, n; + index_t n, dim_max; value_t threshold; + float ratio; coefficient_t modulus; const binomial_coeff_table binomial_coeff; std::vector multiplicative_inverse; mutable std::vector vertices; public: - ripser(compressed_lower_distance_matrix&& _dist, index_t _dim_max, value_t _threshold, + ripser(compressed_lower_distance_matrix&& _dist, index_t _dim_max, value_t _threshold, float _ratio, 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), + : dist(std::move(_dist)), n(dist.size()), + dim_max(std::min(_dim_max, index_t(dist.size() - 2))), threshold(_threshold), + 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 { @@ -624,7 +626,7 @@ public: found_persistence_pair: #ifdef PRINT_PERSISTENCE_PAIRS value_t death = get_diameter(pivot); - if (diameter != death) { + if (death > diameter * ratio) { #ifdef INDICATE_PROGRESS std::cout << "\033[K"; #endif @@ -842,6 +844,8 @@ int main(int argc, char** argv) { index_t dim_max = 1; value_t threshold = std::numeric_limits::infinity(); + float ratio = 1; + #ifdef USE_COEFFICIENTS coefficient_t modulus = 2; #else @@ -862,6 +866,11 @@ int main(int argc, char** argv) { size_t next_pos; threshold = std::stof(parameter, &next_pos); if (next_pos != parameter.size()) print_usage_and_exit(-1); + } else if (arg == "--ratio") { + std::string parameter = std::string(argv[++i]); + size_t next_pos; + ratio = std::stof(parameter, &next_pos); + 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") @@ -915,7 +924,7 @@ int main(int argc, char** argv) { std::cout << "value range: [" << min << "," << max << "]" << std::endl; - ripser(std::move(dist), dim_max, threshold, modulus).compute_barcodes(); + ripser(std::move(dist), dim_max, threshold, ratio, modulus).compute_barcodes(); } void ripser::compute_barcodes() { -- cgit v1.2.3 From ff4b0cc4012e1cdc8ebe0ae5cdc613352b9f8bac Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sat, 19 May 2018 14:17:27 -0400 Subject: consolidation of branches --- ripser.cpp | 447 +++++++++++++++++++++++++++++++++---------------------------- 1 file changed, 244 insertions(+), 203 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 8f2d509..9e291ca 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -141,7 +141,8 @@ public: 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(const diameter_index_t& _diameter_index) : diameter_entry_t(_diameter_index, 1) {} + diameter_entry_t(const diameter_index_t& _diameter_index) + : diameter_entry_t(_diameter_index, 1) {} }; const entry_t& get_entry(const diameter_entry_t& p) { return p.second; } @@ -354,8 +355,8 @@ public: } }; -class ripser { - compressed_lower_distance_matrix dist; +template class ripser { + DistanceMatrix dist; index_t n, dim_max; value_t threshold; float ratio; @@ -363,9 +364,10 @@ class ripser { const binomial_coeff_table binomial_coeff; std::vector multiplicative_inverse; mutable std::vector vertices; + mutable std::vector coface_entries; public: - ripser(compressed_lower_distance_matrix&& _dist, index_t _dim_max, value_t _threshold, float _ratio, + ripser(DistanceMatrix&& _dist, index_t _dim_max, value_t _threshold, float _ratio, coefficient_t _modulus) : dist(std::move(_dist)), n(dist.size()), dim_max(std::min(_dim_max, index_t(dist.size() - 2))), threshold(_threshold), @@ -415,83 +417,84 @@ public: return diam; } - class simplex_coboundary_enumerator { - private: - index_t idx_below, idx_above, v, k; - std::vector vertices; - const diameter_entry_t simplex; - 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) { - parent.get_simplex_vertices(get_index(_simplex), _dim, parent.n, vertices.begin()); - } + class simplex_coboundary_enumerator; - bool has_next() { - while ((v != -1) && (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 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; - return diameter_entry_t(coface_diameter, coface_index, coface_coefficient); - } - }; + void assemble_columns_to_reduce(std::vector& simplices, + std::vector& columns_to_reduce, + hash_map& pivot_column_index, index_t dim); - void compute_barcodes(); + void compute_dim_0_pairs(std::vector& edges, + std::vector& columns_to_reduce) { + union_find dset(n); - 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); + edges = get_edges(); - columns_to_reduce.clear(); + std::sort(edges.rbegin(), edges.rend(), + greater_diameter_or_smaller_index()); -#ifdef INDICATE_PROGRESS - std::cout << "\033[K" - << "assembling " << num_simplices << " columns" << std::flush << "\r"; +#ifdef PRINT_PERSISTENCE_PAIRS + std::cout << "persistence intervals in dim 0:" << std::endl; #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"; + 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 INDICATE_PROGRESS - std::cout << "\033[K" - << "sorting " << num_simplices << " columns" << std::flush << "\r"; +#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 + } - std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), - greater_diameter_or_smaller_index()); -#ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + template + diameter_entry_t add_coboundary_and_get_pivot(Iterator column_begin, Iterator column_end, + coefficient_t factor_column_to_add, +#ifdef ASSEMBLE_REDUCTION_MATRIX + Column& working_reduction_column, #endif + 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); + +#ifdef ASSEMBLE_REDUCTION_MATRIX + working_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()) { + 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, @@ -511,7 +514,8 @@ public: std::vector coface_entries; - for (index_t index_column_to_reduce = 0; index_column_to_reduce < columns_to_reduce.size(); ++index_column_to_reduce) { + 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]; #ifdef ASSEMBLE_REDUCTION_MATRIX @@ -529,14 +533,14 @@ public: #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 << ")" << std::flush << "\r"; + << "reducing column " << index_column_to_reduce + 1 << "/" + << columns_to_reduce.size() << " (diameter " << diameter << ")" + << std::flush << "\r"; #endif index_t index_column_to_add = index_column_to_reduce; - - diameter_entry_t pivot; + 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 @@ -552,8 +556,7 @@ public: bool might_be_apparent_pair = (index_column_to_reduce == index_column_to_add); - do { - + while (true) { #ifdef ASSEMBLE_REDUCTION_MATRIX #ifdef USE_COEFFICIENTS auto reduction_column_begin = reduction_matrix.cbegin(index_column_to_add), @@ -571,39 +574,17 @@ public: auto reduction_column_begin = &reduction_matrix[index_column_to_add], reduction_column_end = &reduction_matrix[index_column_to_add] + 1; #else - auto reduction_column_begin = &columns_to_reduce[index_column_to_add], reduction_column_end = &columns_to_reduce[index_column_to_add] + 1; + auto reduction_column_begin = &columns_to_reduce[index_column_to_add], + reduction_column_end = &columns_to_reduce[index_column_to_add] + 1; #endif #endif - for (auto it = reduction_column_begin; it != reduction_column_end; ++it) { - diameter_entry_t simplex = *it; - set_coefficient(simplex, get_coefficient(simplex) * factor_column_to_add % modulus); - + pivot = add_coboundary_and_get_pivot( + reduction_column_begin, reduction_column_end, factor_column_to_add, #ifdef ASSEMBLE_REDUCTION_MATRIX - working_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); - } - - pivot = get_pivot(working_coboundary, modulus); + working_reduction_column, +#endif + working_coboundary, dim, pivot_column_index, might_be_apparent_pair); if (get_index(pivot) != -1) { auto pair = pivot_column_index.find(get_index(pivot)); @@ -611,69 +592,174 @@ public: if (pair != pivot_column_index.end()) { index_column_to_add = pair->second; factor_column_to_add = modulus - get_coefficient(pivot); - continue; - } - } else { + } else { #ifdef PRINT_PERSISTENCE_PAIRS + value_t death = get_diameter(pivot); + if (diameter != death) { #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; -#endif - std::cout << " [" << diameter << ", )" << std::endl << std::flush; + std::cout << "\033[K"; #endif - break; - } - - found_persistence_pair: -#ifdef PRINT_PERSISTENCE_PAIRS - value_t death = get_diameter(pivot); - if (death > diameter * ratio) { -#ifdef INDICATE_PROGRESS - 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), index_column_to_reduce)); + pivot_column_index.insert( + std::make_pair(get_index(pivot), index_column_to_reduce)); #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_matrix (with a single diagonal 1 entry) -// by reduction_column (possibly with a different entry on the diagonal) + // replace current column of reduction_matrix (with a single diagonal 1 + // entry) by reduction_column (possibly with a different entry on the + // diagonal) #ifdef USE_COEFFICIENTS - reduction_matrix.pop_back(); + reduction_matrix.pop_back(); #else - pop_pivot(working_reduction_column, modulus); + pop_pivot(working_reduction_column, modulus); #endif - while (true) { - diameter_entry_t e = pop_pivot(working_reduction_column, modulus); - if (get_index(e) == -1) break; + while (true) { + diameter_entry_t e = pop_pivot(working_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_matrix.push_back(e); - } + reduction_matrix.push_back(e); + } #else #ifdef USE_COEFFICIENTS - reduction_matrix.pop_back(); - reduction_matrix.push_back(diameter_entry_t(column_to_reduce, inverse)); + reduction_matrix.pop_back(); + reduction_matrix.push_back(diameter_entry_t(column_to_reduce, inverse)); #endif #endif - break; - } while (true); + break; + } + } else { +#ifdef PRINT_PERSISTENCE_PAIRS + std::cout << " [" << diameter << ", )" << std::endl << std::flush; +#endif + break; + } + } } #ifdef INDICATE_PROGRESS std::cout << "\033[K"; #endif } + + std::vector get_edges(); + + void compute_barcodes() { + + std::vector simplices, columns_to_reduce; + + compute_dim_0_pairs(simplices, columns_to_reduce); + + 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(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; + 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) { + 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); + idx_above += binomial_coeff(v, k + 1); + --v; + --k; + assert(k != -1); + } + 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; + return diameter_entry_t(coface_diameter, coface_index, coface_coefficient); + } }; +template <> std::vector ripser::get_edges() { + 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)); + } + 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 +} + enum file_format { LOWER_DISTANCE_MATRIX, UPPER_DISTANCE_MATRIX, @@ -908,74 +994,29 @@ int main(int argc, char** argv) { compressed_lower_distance_matrix dist = read_file(filename ? file_stream : std::cin, format); - std::cout << "distance matrix with " << dist.size() << " points" << std::endl; + value_t min = std::numeric_limits::infinity(), + max = -std::numeric_limits::infinity(), max_finite = max; + int num_edges = 0; - value_t min = std::numeric_limits::infinity(), max = -std::numeric_limits::infinity(); - - for (auto d: dist.distances) { - if (d != std::numeric_limits::infinity() ) { - min = std::min(min, d); - max = std::max(max, d); - } else { - threshold = std::min(threshold, std::numeric_limits::max()); - } + value_t enclosing_radius = std::numeric_limits::infinity(); + for (index_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)); + enclosing_radius = std::min(enclosing_radius, r_i); } - - std::cout << "value range: [" << min << "," << max << "]" - << std::endl; - - ripser(std::move(dist), dim_max, threshold, ratio, 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 (threshold == std::numeric_limits::max()) threshold = enclosing_radius; - 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 (auto d : dist.distances) { + min = std::min(min, d); + max = std::max(max, d); + max_finite = d != std::numeric_limits::infinity() ? std::max(max, d) : max_finite; + if (d <= threshold) ++num_edges; } - for (index_t dim = 1; dim <= dim_max; ++dim) { - hash_map pivot_column_index; - pivot_column_index.reserve(columns_to_reduce.size()); + std::cout << "value range: [" << min << "," << max_finite << "]" << std::endl; - 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); - } - } + std::cout << "distance matrix with " << dist.size() << " points" << std::endl; + ripser(std::move(dist), dim_max, threshold, ratio, modulus) + .compute_barcodes(); } -- cgit v1.2.3 From 535c4115f73fcf052d6206f12fce9d714c716c05 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sat, 19 May 2018 14:26:22 -0400 Subject: consolidation of branches --- ripser.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ripser.cpp b/ripser.cpp index 9e291ca..9e2aa90 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -928,7 +928,7 @@ int main(int argc, char** argv) { file_format format = DISTANCE_MATRIX; index_t dim_max = 1; - value_t threshold = std::numeric_limits::infinity(); + value_t threshold = std::numeric_limits::max(); float ratio = 1; -- cgit v1.2.3