From 257abf0674e7cf81e43c674d7b92a69c67ab7f00 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 15 Sep 2016 16:25:21 +0200 Subject: upgraded XCode project --- ripser.xcodeproj/project.pbxproj | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/ripser.xcodeproj/project.pbxproj b/ripser.xcodeproj/project.pbxproj index 04bc3a6..2aa6179 100644 --- a/ripser.xcodeproj/project.pbxproj +++ b/ripser.xcodeproj/project.pbxproj @@ -88,7 +88,7 @@ 551018401BD63CB300990BFF /* Project object */ = { isa = PBXProject; attributes = { - LastUpgradeCheck = 0700; + LastUpgradeCheck = 0800; ORGANIZATIONNAME = "ulrich-bauer.org"; TargetAttributes = { 551018471BD63CB300990BFF = { @@ -138,8 +138,10 @@ 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_OBJC_ROOT_CLASS = YES_ERROR; + CLANG_WARN_SUSPICIOUS_MOVE = YES; CLANG_WARN_UNREACHABLE_CODE = YES; CLANG_WARN__DUPLICATE_METHOD_MATCH = YES; COPY_PHASE_STRIP = NO; @@ -181,8 +183,10 @@ 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_OBJC_ROOT_CLASS = YES_ERROR; + CLANG_WARN_SUSPICIOUS_MOVE = YES; CLANG_WARN_UNREACHABLE_CODE = YES; CLANG_WARN__DUPLICATE_METHOD_MATCH = YES; COPY_PHASE_STRIP = NO; -- cgit v1.2.3 From 09ad67d33adeb7443685214954c7b447a94ebdc2 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 15 Sep 2016 16:25:39 +0200 Subject: added first version of sparse distance matrix --- ripser.cpp | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/ripser.cpp b/ripser.cpp index 5568f32..b2226fc 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -294,6 +294,27 @@ public: size_t size() const { return rows.size(); } }; +class sparse_distance_matrix { +public: + std::vector> neighbors; + + template + sparse_distance_matrix(const DistanceMatrix& mat, value_t threshold) + : neighbors(mat.size()) + { + + for (index_t i = 1; i < size(); ++i) + for (index_t j = 0; j < size(); ++j) + if (mat(i, j) <= threshold) + neighbors[i].push_back(diameter_index_t(mat(i, j), j)); + } + + value_t operator()(const index_t i, const index_t j) const; + + 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) { @@ -468,6 +489,11 @@ void assemble_columns_to_reduce(std::vector& columns_to_reduce index_t n, value_t threshold, const binomial_coeff_table& binomial_coeff) { index_t num_simplices = binomial_coeff(n, dim + 2); + // iterate over all (previous) columns_to_reduce + // find cofaces, additional vertices + // if additional vertex is larger than all current ones: append to columns_to_reduce + + columns_to_reduce.clear(); #ifdef INDICATE_PROGRESS @@ -584,6 +610,9 @@ void compute_pairs(std::vector& columns_to_reduce, hash_map Date: Fri, 16 Sep 2016 14:47:38 +0200 Subject: intermediate commit --- ripser.cpp | 81 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/ripser.cpp b/ripser.cpp index b2226fc..a477ce6 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -315,6 +315,87 @@ public: }; +template +class set_intersection_enumerator { +public: + InputIteratorCollection first, last; + + set_intersection_enumerator (InputIteratorCollection _first, InputIteratorCollection _last) : first(_first), last(_last) {} + +// template +// set_intersection_enumerator (InputCollection _set) : first(_first), last(_last) {} + + bool has_next() { + for (auto &it0 = first[0], end0 = last[0]; it0 != end0; ++it0) { + auto x = *it0; + for (size_t idx = 1; idx < first.size(); ++idx) { + auto &it = first[idx], end = last[idx]; + while (*it < x) if (++it == end) return false; + if (*it > x) goto continue_outer; + } + return true; + continue_outer: ; + } + return false; + } + + // only valid after has_next() has been called + decltype(*first[0]) next() { + return *first[0]++; + } +}; + + +class simplex_coboundary_enumerator_sparse { +private: + index_t idx, modified_idx, v, k; + const binomial_coeff_table& binomial_coeff; + const sparse_distance_matrix& sparse_dist; + std::vector vertices; + + std::vector::const_reverse_iterator> ii, ee; + set_intersection_enumerator v_intersection; + +public: + simplex_coboundary_enumerator_sparse(index_t _idx, index_t _dim, index_t _n, const binomial_coeff_table& _binomial_coeff, sparse_distance_matrix& _sparse_dist) + : idx(_idx), modified_idx(_idx), v(_n - 1), k(_dim + 1), binomial_coeff(_binomial_coeff), sparse_dist(_sparse_dist), v_intersection(ii, ee) + { + get_simplex_vertices(idx, _dim, _n, _binomial_coeff, std::back_inserter(vertices)); + + for (auto v: vertices) { + ii.push_back(sparse_dist.neighbors[v].rbegin()); + ee.push_back(sparse_dist.neighbors[v].rend()); + } + + vertices.push_back(-1); + + } + + bool has_next() { + return v_intersection.has_next(); + } + + std::pair next() { + + auto v = v_intersection.has_next(); + + while (k > 0 && vertices[k - 1] > v) { + vertices[k] = vertices[k - 1]; + --k; + } + + vertices[k] = v; + + //... + auto result = std::make_pair(make_entry(modified_idx + binomial_coeff(v, k + 1), k & 1 ? -1 : 1), v); + + + return result; + } +}; + + + template <> void compressed_distance_matrix::init_rows() { value_t* pointer = &distances[0]; for (index_t i = 1; i < size(); ++i) { -- cgit v1.2.3 From 01f16a612b8d30a1cc7dc000a5e7d79f00238b2d Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Fri, 16 Sep 2016 19:57:59 +0200 Subject: sparse coboundary enumerator (untested) --- ripser.cpp | 72 ++++++++++++++++++++++++++++++++++---------------------------- 1 file changed, 39 insertions(+), 33 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 28c8b3c..2a62170 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -94,30 +94,36 @@ std::vector multiplicative_inverse_vector(const coefficient_t m) return inverse; } +void get_next_vertex(index_t& v, const index_t idx, const index_t k, const binomial_coeff_table& binomial_coeff) { + + if (binomial_coeff(v, k) > idx) { + index_t count = v; + while (count > 0) { + index_t i = v; + index_t step = count >> 1; + i -= step; + if (binomial_coeff(i, k) > idx) { + v = --i; + count -= step + 1; + } else + count = step; + } + } + assert(binomial_coeff(v, k) <= idx); + assert(binomial_coeff(v + 1, k) > idx); +} + + template -OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t n, +OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v, const binomial_coeff_table& binomial_coeff, OutputIterator out) { - --n; + --v; for (index_t k = dim + 1; k > 0; --k) { - if (binomial_coeff(n, k) > idx) { - index_t count = n; - while (count > 0) { - index_t i = n; - index_t step = count >> 1; - i -= step; - if (binomial_coeff(i, k) > idx) { - n = --i; - count -= step + 1; - } else - count = step; - } - } - assert(binomial_coeff(n, k) <= idx); - assert(binomial_coeff(n + 1, k) > idx); + get_next_vertex(v, idx, k, binomial_coeff); - *out++ = n; - idx -= binomial_coeff(n, k); + *out++ = v; + idx -= binomial_coeff(v, k); } return out; @@ -349,7 +355,7 @@ public: class simplex_coboundary_enumerator_sparse { private: - index_t idx, modified_idx, v, k; + index_t idx_below, idx_above, v, k, w; const binomial_coeff_table& binomial_coeff; const sparse_distance_matrix& sparse_dist; std::vector vertices; @@ -359,16 +365,14 @@ private: public: simplex_coboundary_enumerator_sparse(index_t _idx, index_t _dim, index_t _n, const binomial_coeff_table& _binomial_coeff, sparse_distance_matrix& _sparse_dist) - : idx(_idx), modified_idx(_idx), v(_n - 1), k(_dim + 1), binomial_coeff(_binomial_coeff), sparse_dist(_sparse_dist), v_intersection(ii, ee) + : idx_below(_idx), idx_above(0), v(_n - 1), k(_dim + 1), w(_n - 1), binomial_coeff(_binomial_coeff), sparse_dist(_sparse_dist), v_intersection(ii, ee) { - get_simplex_vertices(idx, _dim, _n, _binomial_coeff, std::back_inserter(vertices)); + get_simplex_vertices(idx_below, _dim, _n, _binomial_coeff, std::back_inserter(vertices)); for (auto v: vertices) { ii.push_back(sparse_dist.neighbors[v].rbegin()); ee.push_back(sparse_dist.neighbors[v].rend()); } - - vertices.push_back(-1); } @@ -378,17 +382,19 @@ public: std::pair next() { - auto v = v_intersection.has_next(); - - while (k > 0 && vertices[k - 1] > v) { - vertices[k] = vertices[k - 1]; - --k; - } + index_t v = get_index(v_intersection.next()); + + while (w > v && idx_below > 0) { + get_next_vertex(w, idx_below, k, binomial_coeff); - vertices[k] = v; + idx_below -= binomial_coeff(w, k); + idx_above += binomial_coeff(w, k + 1); - //... - auto result = std::make_pair(make_entry(modified_idx + binomial_coeff(v, k + 1), k & 1 ? -1 : 1), v); + --k; + assert(k != -1); + } + + auto result = std::make_pair(make_entry(idx_above + binomial_coeff(v, k + 1) + idx_below, k & 1 ? -1 : 1), v); return result; -- cgit v1.2.3 From 5d356d18ab38d85c255927419a45ff6736e1595d Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Fri, 16 Sep 2016 21:14:04 +0200 Subject: sparse coboundary enumerator --- ripser.cpp | 50 +++++++++++++++++++++++++++++++++++++------------- 1 file changed, 37 insertions(+), 13 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 2a62170..85e94f6 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -310,9 +310,9 @@ public: : neighbors(mat.size()) { - for (index_t i = 1; i < size(); ++i) + for (index_t i = 0; i < size(); ++i) for (index_t j = 0; j < size(); ++j) - if (mat(i, j) <= threshold) + if (i != j && mat(i, j) <= threshold) neighbors[i].push_back(diameter_index_t(mat(i, j), j)); } @@ -322,23 +322,42 @@ public: }; -template +template +struct second_argument_greater { + bool operator()(const T &lhs, const T &rhs) const + { + return lhs.second > rhs.second; + } +}; + +template +struct second_argument_equal_to { + bool operator()(const T &lhs, const T &rhs) const + { + return lhs.second == rhs.second; + } +}; + +template class set_intersection_enumerator { public: - InputIteratorCollection first, last; - - set_intersection_enumerator (InputIteratorCollection _first, InputIteratorCollection _last) : first(_first), last(_last) {} + InputIteratorCollection &first, &last; + Compare comp; + Equal equal; + + set_intersection_enumerator (InputIteratorCollection& _first, InputIteratorCollection& _last) : first(_first), last(_last) {} // template // set_intersection_enumerator (InputCollection _set) : first(_first), last(_last) {} bool has_next() { - for (auto &it0 = first[0], end0 = last[0]; it0 != end0; ++it0) { + for (auto &it0 = first[0], &end0 = last[0]; it0 != end0; ++it0) { auto x = *it0; for (size_t idx = 1; idx < first.size(); ++idx) { auto &it = first[idx], end = last[idx]; - while (*it < x) if (++it == end) return false; - if (*it > x) goto continue_outer; + while (comp(*it, x)) if (++it == end) return false; + auto y = *it; + if (!equal(y, x)) goto continue_outer; } return true; continue_outer: ; @@ -361,7 +380,7 @@ private: std::vector vertices; std::vector::const_reverse_iterator> ii, ee; - set_intersection_enumerator v_intersection; + set_intersection_enumerator, second_argument_equal_to> v_intersection; public: simplex_coboundary_enumerator_sparse(index_t _idx, index_t _dim, index_t _n, const binomial_coeff_table& _binomial_coeff, sparse_distance_matrix& _sparse_dist) @@ -373,7 +392,8 @@ public: ii.push_back(sparse_dist.neighbors[v].rbegin()); ee.push_back(sparse_dist.neighbors[v].rend()); } - + + get_next_vertex(w, idx_below, k, binomial_coeff); } bool has_next() { @@ -384,13 +404,17 @@ public: index_t v = get_index(v_intersection.next()); + + while (w > v && idx_below > 0) { - get_next_vertex(w, idx_below, k, binomial_coeff); - idx_below -= binomial_coeff(w, k); idx_above += binomial_coeff(w, k + 1); --k; + + get_next_vertex(w, idx_below, k, binomial_coeff); + + assert(k != -1); } -- cgit v1.2.3 From 29b383df930009dfe43d9629bc666d18345a96ac Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Fri, 16 Sep 2016 21:38:57 +0200 Subject: sparse coboundary enumerator --- ripser.cpp | 29 +++++------------------------ 1 file changed, 5 insertions(+), 24 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 85e94f6..820a3d3 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -94,7 +94,7 @@ std::vector multiplicative_inverse_vector(const coefficient_t m) return inverse; } -void get_next_vertex(index_t& v, const index_t idx, const index_t k, const binomial_coeff_table& binomial_coeff) { +index_t get_next_vertex(index_t& v, const index_t idx, const index_t k, const binomial_coeff_table& binomial_coeff) { if (binomial_coeff(v, k) > idx) { index_t count = v; @@ -111,6 +111,8 @@ void get_next_vertex(index_t& v, const index_t idx, const index_t k, const binom } assert(binomial_coeff(v, k) <= idx); assert(binomial_coeff(v + 1, k) > idx); + + return v; } @@ -321,7 +323,6 @@ public: size_t size() const { return neighbors.size(); } }; - template struct second_argument_greater { bool operator()(const T &lhs, const T &rhs) const @@ -347,9 +348,6 @@ public: set_intersection_enumerator (InputIteratorCollection& _first, InputIteratorCollection& _last) : first(_first), last(_last) {} -// template -// set_intersection_enumerator (InputCollection _set) : first(_first), last(_last) {} - bool has_next() { for (auto &it0 = first[0], &end0 = last[0]; it0 != end0; ++it0) { auto x = *it0; @@ -371,7 +369,6 @@ public: } }; - class simplex_coboundary_enumerator_sparse { private: index_t idx_below, idx_above, v, k, w; @@ -392,8 +389,6 @@ public: ii.push_back(sparse_dist.neighbors[v].rbegin()); ee.push_back(sparse_dist.neighbors[v].rend()); } - - get_next_vertex(w, idx_below, k, binomial_coeff); } bool has_next() { @@ -401,32 +396,18 @@ public: } std::pair next() { - index_t v = get_index(v_intersection.next()); - - - while (w > v && idx_below > 0) { + while (k > 0 && get_next_vertex(w, idx_below, k, binomial_coeff) > v) { idx_below -= binomial_coeff(w, k); idx_above += binomial_coeff(w, k + 1); - --k; - - get_next_vertex(w, idx_below, k, binomial_coeff); - - - assert(k != -1); } - auto result = std::make_pair(make_entry(idx_above + binomial_coeff(v, k + 1) + idx_below, k & 1 ? -1 : 1), v); - - - return result; + return std::make_pair(make_entry(idx_above + binomial_coeff(v, k + 1) + idx_below, k & 1 ? -1 : 1), v); } }; - - template <> void compressed_distance_matrix::init_rows() { value_t* pointer = &distances[0]; for (index_t i = 1; i < size(); ++i) { -- cgit v1.2.3 From 5794ba59f5621c6ec70ffa75a938d61c9342d4e3 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Mon, 19 Sep 2016 10:02:02 +0200 Subject: restructured coboundary enumerator --- ripser.cpp | 292 ++++++++++++++++++++++++++++++++++++++----------------------- 1 file changed, 181 insertions(+), 111 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 208ebb0..80b00d3 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -111,11 +111,10 @@ index_t get_next_vertex(index_t& v, const index_t idx, const index_t k, const bi } assert(binomial_coeff(v, k) <= idx); assert(binomial_coeff(v + 1, k) > idx); - + return v; } - template OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v, const binomial_coeff_table& binomial_coeff, OutputIterator out) { @@ -246,14 +245,26 @@ public: } }; -class simplex_coboundary_enumerator { +template 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 DistanceMatrix& dist; const binomial_coeff_table& binomial_coeff; public: - simplex_coboundary_enumerator(index_t _idx, index_t _dim, index_t _n, const binomial_coeff_table& _binomial_coeff) - : idx_below(_idx), idx_above(0), v(_n - 1), k(_dim + 1), binomial_coeff(_binomial_coeff) {} + simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, index_t _n, + const coefficient_t _modulus, const DistanceMatrix& _dist, + const binomial_coeff_table& _binomial_coeff) + : simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(_n - 1), k(_dim + 1), modulus(_modulus), + binomial_coeff(_binomial_coeff), dist(_dist), vertices(_dim + 1) { + get_simplex_vertices(get_index(_simplex), _dim, _n, binomial_coeff, vertices.begin()); + } bool has_next() { while ((v != -1) && (binomial_coeff(v, k) <= idx_below)) { @@ -267,8 +278,16 @@ public: return v != -1; } - std::pair next() { - auto result = std::make_pair(make_entry(idx_above + binomial_coeff(v, k + 1) + idx_below, k & 1 ? -1 : 1), v); + std::pair next() { + + value_t coface_diameter = get_diameter(simplex); + for (index_t w : vertices) { coface_diameter = std::max(coface_diameter, dist(v, w)); } + + coefficient_t coface_coefficient = (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; + + auto result = std::make_pair( + make_diameter_entry(coface_diameter, idx_above + binomial_coeff(v, k + 1) + idx_below, coface_coefficient), + v); --v; return result; } @@ -308,14 +327,11 @@ public: std::vector> neighbors; template - sparse_distance_matrix(const DistanceMatrix& mat, value_t threshold) - : neighbors(mat.size()) - { + sparse_distance_matrix(const DistanceMatrix& mat, value_t threshold) : neighbors(mat.size()) { for (index_t i = 0; i < size(); ++i) for (index_t j = 0; j < size(); ++j) - if (i != j && mat(i, j) <= threshold) - neighbors[i].push_back(diameter_index_t(mat(i, j), j)); + if (i != j && mat(i, j) <= threshold) neighbors[i].push_back(diameter_index_t(mat(i, j), j)); } value_t operator()(const index_t i, const index_t j) const; @@ -323,88 +339,107 @@ public: size_t size() const { return neighbors.size(); } }; -template -struct second_argument_greater { - bool operator()(const T &lhs, const T &rhs) const - { - return lhs.second > rhs.second; - } +template struct second_argument_greater { + bool operator()(const T& lhs, const T& rhs) const { return lhs.second > rhs.second; } }; -template -struct second_argument_equal_to { - bool operator()(const T &lhs, const T &rhs) const - { - return lhs.second == rhs.second; - } +template struct second_argument_equal_to { + bool operator()(const T& lhs, const T& rhs) const { return lhs.second == rhs.second; } }; -template +template class set_intersection_enumerator { public: - InputIteratorCollection &first, &last; + InputIteratorCollection &ii, ⅇ Compare comp; Equal equal; - - set_intersection_enumerator (InputIteratorCollection& _first, InputIteratorCollection& _last) : first(_first), last(_last) {} - - bool has_next() { - for (auto &it0 = first[0], &end0 = last[0]; it0 != end0; ++it0) { - auto x = *it0; - for (size_t idx = 1; idx < first.size(); ++idx) { - auto &it = first[idx], end = last[idx]; - while (comp(*it, x)) if (++it == end) return false; + Representative rep; + + set_intersection_enumerator(InputIteratorCollection& _first, InputIteratorCollection& _last) + : ii(_first), ee(_last) {} + + bool has_next() { + for (auto &it0 = ii[0], &end0 = ee[0]; it0 != end0; ++it0) { + auto x = *it0; + for (size_t idx = 1; idx < ii.size(); ++idx) { + auto &it = ii[idx], end = ee[idx]; + while (comp(*it, x)) + if (++it == end) return false; auto y = *it; - if (!equal(y, x)) goto continue_outer; - } - return true; - continue_outer: ; - } - return false; - } - - // only valid after has_next() has been called - decltype(*first[0]) next() { - return *first[0]++; - } + if (!equal(y, x)) + goto continue_outer; + else + x = rep(x, y); + } + return true; + continue_outer:; + } + return false; + } + + // only valid after has_next() has been called + decltype(*ii[0]) next() { return *ii[0]++; } }; class simplex_coboundary_enumerator_sparse { private: - index_t idx_below, idx_above, v, k, w; + index_t idx_below, idx_above, v, k, max_vertex_below; const binomial_coeff_table& binomial_coeff; - const sparse_distance_matrix& sparse_dist; - std::vector vertices; - - std::vector::const_reverse_iterator> ii, ee; - set_intersection_enumerator, second_argument_equal_to> v_intersection; + const sparse_distance_matrix& sparse_dist; + std::vector vertices; + + const value_t simplex_diameter; + + std::vector::const_reverse_iterator> ii, ee; + + diameter_index_t x; public: - simplex_coboundary_enumerator_sparse(index_t _idx, index_t _dim, index_t _n, const binomial_coeff_table& _binomial_coeff, sparse_distance_matrix& _sparse_dist) - : idx_below(_idx), idx_above(0), v(_n - 1), k(_dim + 1), w(_n - 1), binomial_coeff(_binomial_coeff), sparse_dist(_sparse_dist), v_intersection(ii, ee) - { - get_simplex_vertices(idx_below, _dim, _n, _binomial_coeff, std::back_inserter(vertices)); - - for (auto v: vertices) { - ii.push_back(sparse_dist.neighbors[v].rbegin()); - ee.push_back(sparse_dist.neighbors[v].rend()); - } - } + simplex_coboundary_enumerator_sparse(diameter_index_t _simplex, index_t _dim, index_t _n, + const binomial_coeff_table& _binomial_coeff, + const sparse_distance_matrix& _sparse_dist) + : simplex_diameter(get_diameter(_simplex)), idx_below(get_index(_simplex)), idx_above(0), v(_n - 1), + k(_dim + 1), max_vertex_below(_n - 1), binomial_coeff(_binomial_coeff), sparse_dist(_sparse_dist) { + get_simplex_vertices(idx_below, _dim, _n, _binomial_coeff, std::back_inserter(vertices)); + + for (auto v : vertices) { + ii.push_back(sparse_dist.neighbors[v].rbegin()); + ee.push_back(sparse_dist.neighbors[v].rend()); + } + } bool has_next() { - return v_intersection.has_next(); + for (auto &it0 = ii[0], &end0 = ee[0]; it0 != end0; ++it0) { + x = *it0; + for (size_t idx = 1; idx < ii.size(); ++idx) { + auto &it = ii[idx], end = ee[idx]; + while (get_index(*it) > get_index(x)) + if (++it == end) return false; + auto y = *it; + if (get_index(y) != get_index(x)) + goto continue_outer; + else + x = std::min(x, y, greater_diameter_or_smaller_index()); + } + return true; + continue_outer:; + } + return false; } - std::pair next() { - index_t v = get_index(v_intersection.next()); - - while (k > 0 && get_next_vertex(w, idx_below, k, binomial_coeff) > v) { - idx_below -= binomial_coeff(w, k); - idx_above += binomial_coeff(w, k + 1); + std::pair next() { + index_t covertex = get_index(x); + + while (k > 0 && get_next_vertex(max_vertex_below, idx_below, k, binomial_coeff) > covertex) { + idx_below -= binomial_coeff(max_vertex_below, k); + idx_above += binomial_coeff(max_vertex_below, k + 1); --k; } - - return std::make_pair(make_entry(idx_above + binomial_coeff(v, k + 1) + idx_below, k & 1 ? -1 : 1), v); + + return std::make_pair(make_diameter_entry(get_diameter(x), + idx_above + binomial_coeff(covertex, k + 1) + idx_below, + k & 1 ? -1 : 1), + covertex); } }; @@ -582,10 +617,60 @@ void assemble_columns_to_reduce(std::vector& columns_to_reduce index_t n, value_t threshold, const binomial_coeff_table& binomial_coeff) { index_t num_simplices = binomial_coeff(n, dim + 2); - // iterate over all (previous) columns_to_reduce - // find cofaces, additional vertices - // if additional vertex is larger than all current ones: append to columns_to_reduce + 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 = comp.diameter(index); + if (diameter <= threshold) columns_to_reduce.push_back(std::make_pair(diameter, index)); + } + } + +#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 assemble_columns_to_reduce_sparse(std::vector& columns_to_reduce, + hash_map& pivot_column_index, const Comparator& comp, + index_t dim, index_t n, value_t threshold, + const binomial_coeff_table& binomial_coeff, + const sparse_distance_matrix& sparse_dist) { + index_t num_simplices = binomial_coeff(n, dim + 2); + + // iterate over all (previous) columns_to_reduce + // find cofaces, additional vertices + // if additional vertex is larger than all current ones: append to columns_to_reduce + + std::vector previous_columns_to_reduce; + previous_columns_to_reduce.swap(columns_to_reduce); + for (diameter_index_t simplex : previous_columns_to_reduce) { + // coface_entries.clear(); + simplex_coboundary_enumerator_sparse cofaces(simplex, dim, n, binomial_coeff, sparse_dist); + + while (cofaces.has_next()) { + auto coface_descriptor = cofaces.next(); + auto coface = coface_descriptor.first; + index_t covertex = get_index(coface_descriptor.second); + index_t coface_index = get_index(coface); + value_t coface_diameter = get_diameter(coface_descriptor.first); + // for (index_t v : vertices) { coface_diameter = std::max(coface_diameter, dist(v, covertex)); } + assert(comp.diameter(coface_index) == coface_diameter); + } + } columns_to_reduce.clear(); @@ -615,9 +700,9 @@ void assemble_columns_to_reduce(std::vector& columns_to_reduce template void compute_pairs(std::vector& columns_to_reduce, hash_map& pivot_column_index, - const DistanceMatrix& dist, const ComparatorCofaces& comp, const Comparator& comp_prev, index_t dim, - index_t n, value_t threshold, coefficient_t modulus, - const std::vector& multiplicative_inverse, + index_t dim, index_t n, value_t threshold, coefficient_t modulus, + const std::vector& multiplicative_inverse, const DistanceMatrix& dist, + const ComparatorCofaces& comp, const Comparator& comp_prev, const binomial_coeff_table& binomial_coeff) { #ifdef PRINT_PERSISTENCE_PAIRS @@ -682,51 +767,36 @@ void compute_pairs(std::vector& columns_to_reduce, hash_map cofaces(simplex, dim, n, modulus, dist, binomial_coeff); + while (cofaces.has_next()) { - auto coface_descriptor = cofaces.next(); - entry_t coface = coface_descriptor.first; - index_t covertex = coface_descriptor.second; - index_t coface_index = get_index(coface); - value_t coface_diameter = get_diameter(simplex); - for (index_t v : vertices) { coface_diameter = std::max(coface_diameter, dist(v, covertex)); } - assert(comp.diameter(coface_index) == coface_diameter); - - if (coface_diameter <= threshold) { - coefficient_t coface_coefficient = - (get_coefficient(coface) + modulus) * simplex_coefficient % modulus; - assert(coface_coefficient >= 0); - - diameter_entry_t coface_entry = - make_diameter_entry(coface_diameter, coface_index, coface_coefficient); - coface_entries.push_back(coface_entry); - - if (might_be_apparent_pair && (get_diameter(simplex) == coface_diameter)) { - if (pivot_column_index.find(coface_index) == pivot_column_index.end()) { - pivot = coface_entry; + diameter_entry_t coface = cofaces.next().first; + assert(comp.diameter(get_index(coface)) == get_diameter(coface)); + + 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; @@ -1064,8 +1134,8 @@ int main(int argc, char** argv) { hash_map pivot_column_index; pivot_column_index.reserve(columns_to_reduce.size()); - compute_pairs(columns_to_reduce, pivot_column_index, dist, comp, comp_prev, dim, n, threshold, modulus, - multiplicative_inverse, binomial_coeff); + compute_pairs(columns_to_reduce, pivot_column_index, dim, n, threshold, modulus, multiplicative_inverse, dist, + comp, comp_prev, binomial_coeff); if (dim < dim_max) { assemble_columns_to_reduce(columns_to_reduce, pivot_column_index, comp, dim, n, threshold, binomial_coeff); -- cgit v1.2.3 From 51b07afc8e3d68d2c8dac8ef2916672ecffd9bfa Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Mon, 19 Sep 2016 18:05:09 +0200 Subject: cleanup --- ripser.cpp | 112 ++++++++++++++++++++++++++++++------------------------------- 1 file changed, 55 insertions(+), 57 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 80b00d3..c2213df 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -178,7 +178,11 @@ template struct smaller_index { bool operator()(const Entry& a, const Entry& b) { return get_index(a) < get_index(b); } }; -typedef std::pair diameter_index_t; +class diameter_index_t: public std::pair { +public: + diameter_index_t() : std::pair() {} + diameter_index_t(std::pair p) : std::pair(p) {} +}; value_t get_diameter(diameter_index_t i) { return i.first; } index_t get_index(diameter_index_t i) { return i.second; } @@ -187,6 +191,12 @@ public: diameter_entry_t(std::pair p) : std::pair(p) {} diameter_entry_t(entry_t e) : std::pair(0, e) {} diameter_entry_t() : diameter_entry_t(0) {} + diameter_entry_t(value_t _diameter, index_t _index, coefficient_t _coefficient) + : std::pair(_diameter, make_entry(_index, _coefficient)) {} + diameter_entry_t(diameter_index_t _diameter_index, coefficient_t _coefficient) + : std::pair(get_diameter(_diameter_index), + make_entry(get_index(_diameter_index), _coefficient)) {} + diameter_entry_t(diameter_index_t _diameter_index) : diameter_entry_t(_diameter_index, 1) {} }; const entry_t& get_entry(const diameter_entry_t& p) { return p.second; } @@ -195,12 +205,6 @@ const index_t get_index(const diameter_entry_t& p) { return get_index(get_entry( 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); } -diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index, coefficient_t _coefficient) { - return std::make_pair(_diameter, make_entry(_index, _coefficient)); -} -diameter_entry_t make_diameter_entry(diameter_index_t _diameter_index, coefficient_t _coefficient) { - return std::make_pair(get_diameter(_diameter_index), make_entry(get_index(_diameter_index), _coefficient)); -} template struct greater_diameter_or_smaller_index { bool operator()(const Entry& a, const Entry& b) { @@ -278,18 +282,16 @@ public: return v != -1; } - std::pair next() { + index_t next_index() { return idx_above + binomial_coeff(v--, k + 1) + idx_below; } + + 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)); } coefficient_t coface_coefficient = (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; - auto result = std::make_pair( - make_diameter_entry(coface_diameter, idx_above + binomial_coeff(v, k + 1) + idx_below, coface_coefficient), - v); - --v; - return result; + return diameter_entry_t(coface_diameter, idx_above + binomial_coeff(v, k + 1) + idx_below, coface_coefficient); } }; @@ -388,18 +390,20 @@ private: const sparse_distance_matrix& sparse_dist; std::vector vertices; - const value_t simplex_diameter; + const diameter_entry_t simplex; std::vector::const_reverse_iterator> ii, ee; diameter_index_t x; + const coefficient_t modulus; + public: - simplex_coboundary_enumerator_sparse(diameter_index_t _simplex, index_t _dim, index_t _n, - const binomial_coeff_table& _binomial_coeff, - const sparse_distance_matrix& _sparse_dist) - : simplex_diameter(get_diameter(_simplex)), idx_below(get_index(_simplex)), idx_above(0), v(_n - 1), - k(_dim + 1), max_vertex_below(_n - 1), binomial_coeff(_binomial_coeff), sparse_dist(_sparse_dist) { + simplex_coboundary_enumerator_sparse(const diameter_entry_t _simplex, index_t _dim, index_t _n, + const coefficient_t _modulus, const sparse_distance_matrix& _sparse_dist, + const binomial_coeff_table& _binomial_coeff) + : simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(_n - 1), k(_dim + 1), + max_vertex_below(_n - 1), modulus(_modulus), sparse_dist(_sparse_dist), binomial_coeff(_binomial_coeff) { get_simplex_vertices(idx_below, _dim, _n, _binomial_coeff, std::back_inserter(vertices)); for (auto v : vertices) { @@ -427,19 +431,21 @@ public: return false; } - std::pair next() { - index_t covertex = get_index(x); + diameter_entry_t next() { - while (k > 0 && get_next_vertex(max_vertex_below, idx_below, k, binomial_coeff) > covertex) { + while (k > 0 && get_next_vertex(max_vertex_below, idx_below, k, binomial_coeff) > get_index(x)) { idx_below -= binomial_coeff(max_vertex_below, k); idx_above += binomial_coeff(max_vertex_below, k + 1); --k; } - return std::make_pair(make_diameter_entry(get_diameter(x), - idx_above + binomial_coeff(covertex, k + 1) + idx_below, - k & 1 ? -1 : 1), - covertex); + value_t coface_diameter = std::max(get_diameter(simplex), get_diameter(x)); + + coefficient_t coface_coefficient = (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; + + return diameter_entry_t(get_diameter(x), + idx_above + binomial_coeff(get_index(x), k + 1) + idx_below, + k & 1 ? -1 : 1); } }; @@ -646,10 +652,9 @@ void assemble_columns_to_reduce(std::vector& columns_to_reduce template void assemble_columns_to_reduce_sparse(std::vector& columns_to_reduce, hash_map& pivot_column_index, const Comparator& comp, - index_t dim, index_t n, value_t threshold, - const binomial_coeff_table& binomial_coeff, - const sparse_distance_matrix& sparse_dist) { - index_t num_simplices = binomial_coeff(n, dim + 2); + index_t dim, index_t n, value_t threshold, const coefficient_t modulus, + const sparse_distance_matrix& sparse_dist, + const binomial_coeff_table& binomial_coeff) { // iterate over all (previous) columns_to_reduce // find cofaces, additional vertices @@ -659,16 +664,12 @@ void assemble_columns_to_reduce_sparse(std::vector& columns_to previous_columns_to_reduce.swap(columns_to_reduce); for (diameter_index_t simplex : previous_columns_to_reduce) { // coface_entries.clear(); - simplex_coboundary_enumerator_sparse cofaces(simplex, dim, n, binomial_coeff, sparse_dist); + simplex_coboundary_enumerator_sparse cofaces(simplex, dim, n, modulus, sparse_dist, binomial_coeff); while (cofaces.has_next()) { - auto coface_descriptor = cofaces.next(); - auto coface = coface_descriptor.first; - index_t covertex = get_index(coface_descriptor.second); - index_t coface_index = get_index(coface); - value_t coface_diameter = get_diameter(coface_descriptor.first); - // for (index_t v : vertices) { coface_diameter = std::max(coface_diameter, dist(v, covertex)); } - assert(comp.diameter(coface_index) == coface_diameter); + auto coface = cofaces.next(); + + columns_to_reduce.push_back(std::make_pair(get_diameter(coface), get_index(coface))); } } @@ -679,12 +680,12 @@ void assemble_columns_to_reduce_sparse(std::vector& columns_to << "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 = comp.diameter(index); - 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 = comp.diameter(index); +// if (diameter <= threshold) columns_to_reduce.push_back(std::make_pair(diameter, index)); +// } +// } #ifdef INDICATE_PROGRESS std::cout << "\033[K" @@ -745,15 +746,15 @@ void compute_pairs(std::vector& columns_to_reduce, hash_map& columns_to_reduce, hash_map cofaces(simplex, dim, n, modulus, dist, binomial_coeff); while (cofaces.has_next()) { - diameter_entry_t coface = cofaces.next().first; + diameter_entry_t coface = cofaces.next(); assert(comp.diameter(get_index(coface)) == get_diameter(coface)); if (get_diameter(coface) <= threshold) { @@ -848,20 +849,17 @@ void compute_pairs(std::vector& columns_to_reduce, hash_map 0); -#else - const coefficient_t coefficient = 1; + set_coefficient(e, inverse * get_coefficient(e) % modulus); + assert(get_coefficient(e) > 0); #endif - reduction_matrix.push_back(make_diameter_entry(get_diameter(e), index, coefficient)); + reduction_matrix.push_back(e); } #else #ifdef USE_COEFFICIENTS reduction_coefficients.pop_back(); - reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, inverse)); + reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, inverse)); #endif #endif break; @@ -1097,7 +1095,7 @@ int main(int argc, char** argv) { rips_filtration_comparator comp(dist, 1, binomial_coeff); for (index_t index = binomial_coeff(n, 2); index-- > 0;) { value_t diameter = comp.diameter(index); - if (diameter <= threshold) edges.push_back(diameter_index_t(diameter, index)); + if (diameter <= threshold) edges.push_back(std::make_pair(diameter, index)); } std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index()); -- cgit v1.2.3 From ad46246d05aaef70d523b96d483e5a6407404cb1 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Mon, 19 Sep 2016 18:22:37 +0200 Subject: small bugs fixed --- ripser.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index c2213df..1405865 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -291,7 +291,7 @@ public: coefficient_t coface_coefficient = (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; - return diameter_entry_t(coface_diameter, idx_above + binomial_coeff(v, k + 1) + idx_below, coface_coefficient); + return diameter_entry_t(coface_diameter, idx_above + binomial_coeff(v--, k + 1) + idx_below, coface_coefficient); } }; @@ -443,9 +443,9 @@ public: coefficient_t coface_coefficient = (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; - return diameter_entry_t(get_diameter(x), + return diameter_entry_t(coface_diameter, idx_above + binomial_coeff(get_index(x), k + 1) + idx_below, - k & 1 ? -1 : 1); + coface_coefficient); } }; -- cgit v1.2.3 From 42a60a47c9bf12c5249fdea9b354fe352170950f Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Tue, 20 Sep 2016 09:35:45 +0200 Subject: bug fix in sparse coboundary enumerator --- ripser.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ripser.cpp b/ripser.cpp index 1405865..e16679c 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -333,7 +333,7 @@ public: for (index_t i = 0; i < size(); ++i) for (index_t j = 0; j < size(); ++j) - if (i != j && mat(i, j) <= threshold) neighbors[i].push_back(diameter_index_t(mat(i, j), j)); + if (i != j && mat(i, j) <= threshold) neighbors[i].push_back(std::make_pair(mat(i, j), j)); } value_t operator()(const index_t i, const index_t j) const; @@ -432,6 +432,7 @@ public: } diameter_entry_t next() { + ++ii[0]; while (k > 0 && get_next_vertex(max_vertex_below, idx_below, k, binomial_coeff) > get_index(x)) { idx_below -= binomial_coeff(max_vertex_below, k); -- cgit v1.2.3 From ab9fadd86b09a34e159fa12bdf0ab285ca252c3f Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Tue, 20 Sep 2016 11:08:59 +0200 Subject: fixes to sparse cobounrady computation --- ripser.cpp | 44 +++++++++++++++++++++++++--------------- ripser.xcodeproj/project.pbxproj | 11 +++++----- 2 files changed, 34 insertions(+), 21 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index e16679c..5bbfe5c 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -168,7 +168,7 @@ typedef index_t entry_t; const index_t get_index(entry_t i) { return i; } index_t get_coefficient(entry_t i) { return 1; } entry_t make_entry(index_t _index, coefficient_t _value) { return entry_t(_index); } -void set_coefficient(index_t& e, const coefficient_t c) { e = c; } +void set_coefficient(index_t& e, const coefficient_t c) { } #endif @@ -336,8 +336,6 @@ public: if (i != j && mat(i, j) <= threshold) neighbors[i].push_back(std::make_pair(mat(i, j), j)); } - value_t operator()(const index_t i, const index_t j) const; - size_t size() const { return neighbors.size(); } }; @@ -383,7 +381,8 @@ public: decltype(*ii[0]) next() { return *ii[0]++; } }; -class simplex_coboundary_enumerator_sparse { +template<> +class simplex_coboundary_enumerator { private: index_t idx_below, idx_above, v, k, max_vertex_below; const binomial_coeff_table& binomial_coeff; @@ -399,7 +398,7 @@ private: const coefficient_t modulus; public: - simplex_coboundary_enumerator_sparse(const diameter_entry_t _simplex, index_t _dim, index_t _n, + simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, index_t _n, const coefficient_t _modulus, const sparse_distance_matrix& _sparse_dist, const binomial_coeff_table& _binomial_coeff) : simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(_n - 1), k(_dim + 1), @@ -661,11 +660,19 @@ void assemble_columns_to_reduce_sparse(std::vector& columns_to // find cofaces, additional vertices // if additional vertex is larger than all current ones: append to columns_to_reduce +#ifdef INDICATE_PROGRESS + std::cout << "\033[K" + << "assembling columns" << std::flush << "\r"; +#endif + std::vector previous_columns_to_reduce; previous_columns_to_reduce.swap(columns_to_reduce); + + columns_to_reduce.clear(); + for (diameter_index_t simplex : previous_columns_to_reduce) { // coface_entries.clear(); - simplex_coboundary_enumerator_sparse cofaces(simplex, dim, n, modulus, sparse_dist, binomial_coeff); + simplex_coboundary_enumerator cofaces(simplex, dim, n, modulus, sparse_dist, binomial_coeff); while (cofaces.has_next()) { auto coface = cofaces.next(); @@ -674,12 +681,6 @@ void assemble_columns_to_reduce_sparse(std::vector& columns_to } } - 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()) { @@ -690,7 +691,7 @@ void assemble_columns_to_reduce_sparse(std::vector& columns_to #ifdef INDICATE_PROGRESS std::cout << "\033[K" - << "sorting " << num_simplices << " columns" << std::flush << "\r"; + << "sorting " << columns_to_reduce.size() << " columns" << std::flush << "\r"; #endif std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), @@ -700,10 +701,12 @@ void assemble_columns_to_reduce_sparse(std::vector& columns_to #endif } -template +template void compute_pairs(std::vector& columns_to_reduce, hash_map& pivot_column_index, index_t dim, index_t n, value_t threshold, coefficient_t modulus, const std::vector& multiplicative_inverse, const DistanceMatrix& dist, +// const FullDistanceMatrix& dist_full, const ComparatorCofaces& comp, const Comparator& comp_prev, const binomial_coeff_table& binomial_coeff) { @@ -783,15 +786,22 @@ void compute_pairs(std::vector& columns_to_reduce, hash_map cofaces(simplex, dim, n, modulus, dist, binomial_coeff); +// simplex_coboundary_enumerator cofaces_full(simplex, dim, n, modulus, dist_full, binomial_coeff); while (cofaces.has_next()) { +// assert(cofaces_full.has_next()); diameter_entry_t coface = cofaces.next(); - assert(comp.diameter(get_index(coface)) == get_diameter(coface)); +// , coface_full = cofaces_full.next(); +// value_t diameter = comp.diameter(get_index(coface)); +// assert(comp.diameter(get_index(coface)) == get_diameter(coface)); if (get_diameter(coface) <= threshold) { coface_entries.push_back(coface); @@ -1083,6 +1093,8 @@ int main(int argc, char** argv) { auto value_range = std::minmax_element(dist.distances.begin(), dist.distances.end()); std::cout << "value range: [" << *value_range.first << "," << *value_range.second << "]" << std::endl; + sparse_distance_matrix sparse_dist(dist, threshold); + dim_max = std::min(dim_max, n - 2); binomial_coeff_table binomial_coeff(n, dim_max + 2); @@ -1133,7 +1145,7 @@ int main(int argc, char** argv) { hash_map pivot_column_index; pivot_column_index.reserve(columns_to_reduce.size()); - compute_pairs(columns_to_reduce, pivot_column_index, dim, n, threshold, modulus, multiplicative_inverse, dist, + compute_pairs(columns_to_reduce, pivot_column_index, dim, n, threshold, modulus, multiplicative_inverse, sparse_dist, //dist, comp, comp_prev, binomial_coeff); if (dim < dim_max) { diff --git a/ripser.xcodeproj/project.pbxproj b/ripser.xcodeproj/project.pbxproj index 2aa6179..24c7d8b 100644 --- a/ripser.xcodeproj/project.pbxproj +++ b/ripser.xcodeproj/project.pbxproj @@ -129,7 +129,7 @@ isa = XCBuildConfiguration; buildSettings = { ALWAYS_SEARCH_USER_PATHS = NO; - CLANG_CXX_LANGUAGE_STANDARD = "gnu++0x"; + CLANG_CXX_LANGUAGE_STANDARD = "c++0x"; CLANG_CXX_LIBRARY = "libc++"; CLANG_ENABLE_MODULES = YES; CLANG_ENABLE_OBJC_ARC = YES; @@ -148,7 +148,7 @@ DEBUG_INFORMATION_FORMAT = dwarf; ENABLE_STRICT_OBJC_MSGSEND = YES; ENABLE_TESTABILITY = YES; - GCC_C_LANGUAGE_STANDARD = gnu99; + GCC_C_LANGUAGE_STANDARD = c11; GCC_DYNAMIC_NO_PIC = NO; GCC_NO_COMMON_BLOCKS = YES; GCC_OPTIMIZATION_LEVEL = 0; @@ -174,7 +174,7 @@ isa = XCBuildConfiguration; buildSettings = { ALWAYS_SEARCH_USER_PATHS = NO; - CLANG_CXX_LANGUAGE_STANDARD = "gnu++0x"; + CLANG_CXX_LANGUAGE_STANDARD = "c++0x"; CLANG_CXX_LIBRARY = "libc++"; CLANG_ENABLE_MODULES = YES; CLANG_ENABLE_OBJC_ARC = YES; @@ -193,7 +193,7 @@ DEBUG_INFORMATION_FORMAT = "dwarf-with-dsym"; ENABLE_NS_ASSERTIONS = NO; ENABLE_STRICT_OBJC_MSGSEND = YES; - GCC_C_LANGUAGE_STANDARD = gnu99; + GCC_C_LANGUAGE_STANDARD = c11; GCC_NO_COMMON_BLOCKS = YES; GCC_WARN_64_TO_32_BIT_CONVERSION = YES; GCC_WARN_ABOUT_RETURN_TYPE = YES_ERROR; @@ -216,7 +216,8 @@ _FILE_FORMAT_DISTANCE_MATRIX, _FILE_FORMAT_LOWER_DISTANCE_MATRIX, _FILE_FORMAT_POINT_CLOUD, - USE_COEFFICIENTS, + _USE_COEFFICIENTS, + INDICATE_PROGRESS, ); PRODUCT_NAME = "$(TARGET_NAME)"; }; -- cgit v1.2.3 From c6d9d81dc7f445a0bcbcc7e943b514edad448389 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Tue, 20 Sep 2016 13:36:51 +0200 Subject: finished first version of sparse matrix support --- ripser.cpp | 96 +++++++++++++++++++++++--------------------------------------- 1 file changed, 35 insertions(+), 61 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 5bbfe5c..9f4e5e5 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -347,39 +347,6 @@ template struct second_argument_equal_to { bool operator()(const T& lhs, const T& rhs) const { return lhs.second == rhs.second; } }; -template -class set_intersection_enumerator { -public: - InputIteratorCollection &ii, ⅇ - Compare comp; - Equal equal; - Representative rep; - - set_intersection_enumerator(InputIteratorCollection& _first, InputIteratorCollection& _last) - : ii(_first), ee(_last) {} - - bool has_next() { - for (auto &it0 = ii[0], &end0 = ee[0]; it0 != end0; ++it0) { - auto x = *it0; - for (size_t idx = 1; idx < ii.size(); ++idx) { - auto &it = ii[idx], end = ee[idx]; - while (comp(*it, x)) - if (++it == end) return false; - auto y = *it; - if (!equal(y, x)) - goto continue_outer; - else - x = rep(x, y); - } - return true; - continue_outer:; - } - return false; - } - - // only valid after has_next() has been called - decltype(*ii[0]) next() { return *ii[0]++; } -}; template<> class simplex_coboundary_enumerator { @@ -411,7 +378,7 @@ public: } } - bool has_next() { + bool has_next(bool all_cofaces = true) { for (auto &it0 = ii[0], &end0 = ee[0]; it0 != end0; ++it0) { x = *it0; for (size_t idx = 1; idx < ii.size(); ++idx) { @@ -424,7 +391,7 @@ public: else x = std::min(x, y, greater_diameter_or_smaller_index()); } - return true; + return all_cofaces || !(k > 0 && get_next_vertex(max_vertex_below, idx_below, k, binomial_coeff) > get_index(x)); continue_outer:; } return false; @@ -634,6 +601,11 @@ void assemble_columns_to_reduce(std::vector& columns_to_reduce if (pivot_column_index.find(index) == pivot_column_index.end()) { value_t diameter = comp.diameter(index); 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"; +#endif } } @@ -649,11 +621,15 @@ void assemble_columns_to_reduce(std::vector& columns_to_reduce #endif } + + template -void assemble_columns_to_reduce_sparse(std::vector& columns_to_reduce, - hash_map& pivot_column_index, const Comparator& comp, +void assemble_columns_to_reduce(std::vector& simplices, + std::vector& columns_to_reduce, + hash_map& pivot_column_index, + const sparse_distance_matrix& sparse_dist, + const Comparator& comp, index_t dim, index_t n, value_t threshold, const coefficient_t modulus, - const sparse_distance_matrix& sparse_dist, const binomial_coeff_table& binomial_coeff) { // iterate over all (previous) columns_to_reduce @@ -665,30 +641,26 @@ void assemble_columns_to_reduce_sparse(std::vector& columns_to << "assembling columns" << std::flush << "\r"; #endif - std::vector previous_columns_to_reduce; - previous_columns_to_reduce.swap(columns_to_reduce); - columns_to_reduce.clear(); + + std::vector next_simplices; - for (diameter_index_t simplex : previous_columns_to_reduce) { + for (diameter_index_t simplex : simplices) { // coface_entries.clear(); - simplex_coboundary_enumerator cofaces(simplex, dim, n, modulus, sparse_dist, binomial_coeff); + simplex_coboundary_enumerator cofaces(simplex, dim, n, modulus, sparse_dist, binomial_coeff); - while (cofaces.has_next()) { + while (cofaces.has_next(false)) { auto coface = cofaces.next(); + + next_simplices.push_back(std::make_pair(get_diameter(coface), get_index(coface))); - columns_to_reduce.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))); } } - - -// for (index_t index = 0; index < num_simplices; ++index) { -// if (pivot_column_index.find(index) == pivot_column_index.end()) { -// value_t diameter = comp.diameter(index); -// if (diameter <= threshold) columns_to_reduce.push_back(std::make_pair(diameter, index)); -// } -// } - + + simplices.swap(next_simplices); + #ifdef INDICATE_PROGRESS std::cout << "\033[K" << "sorting " << columns_to_reduce.size() << " columns" << std::flush << "\r"; @@ -1102,14 +1074,16 @@ int main(int argc, char** argv) { std::vector columns_to_reduce; + std::vector simplices , &edges = simplices; + rips_filtration_comparator comp(dist, 1, binomial_coeff); + for (index_t index = binomial_coeff(n, 2); index-- > 0;) { + value_t diameter = comp.diameter(index); + if (diameter <= threshold) edges.push_back(std::make_pair(diameter, index)); + } + { union_find dset(n); - std::vector edges; - rips_filtration_comparator comp(dist, 1, binomial_coeff); - for (index_t index = binomial_coeff(n, 2); index-- > 0;) { - value_t diameter = comp.diameter(index); - 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 @@ -1149,7 +1123,7 @@ int main(int argc, char** argv) { comp, comp_prev, binomial_coeff); if (dim < dim_max) { - assemble_columns_to_reduce(columns_to_reduce, pivot_column_index, comp, dim, n, threshold, binomial_coeff); + assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, sparse_dist, comp, dim, n, threshold, modulus, binomial_coeff); } } } -- cgit v1.2.3 From f31f1d94bc2367c27e306562c2f0492d2f0e87ab Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 22 Sep 2016 13:53:04 +0200 Subject: cleaned up comments --- ripser.cpp | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 9f4e5e5..713c70b 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -632,10 +632,6 @@ void assemble_columns_to_reduce(std::vector& simplices, index_t dim, index_t n, value_t threshold, const coefficient_t modulus, const binomial_coeff_table& binomial_coeff) { - // iterate over all (previous) columns_to_reduce - // find cofaces, additional vertices - // if additional vertex is larger than all current ones: append to columns_to_reduce - #ifdef INDICATE_PROGRESS std::cout << "\033[K" << "assembling columns" << std::flush << "\r"; @@ -646,7 +642,6 @@ void assemble_columns_to_reduce(std::vector& simplices, std::vector next_simplices; for (diameter_index_t simplex : simplices) { - // coface_entries.clear(); simplex_coboundary_enumerator cofaces(simplex, dim, n, modulus, sparse_dist, binomial_coeff); while (cofaces.has_next(false)) { @@ -673,12 +668,11 @@ void assemble_columns_to_reduce(std::vector& simplices, #endif } -template void compute_pairs(std::vector& columns_to_reduce, hash_map& pivot_column_index, index_t dim, index_t n, value_t threshold, coefficient_t modulus, const std::vector& multiplicative_inverse, const DistanceMatrix& dist, -// const FullDistanceMatrix& dist_full, const ComparatorCofaces& comp, const Comparator& comp_prev, const binomial_coeff_table& binomial_coeff) { @@ -758,22 +752,14 @@ void compute_pairs(std::vector& columns_to_reduce, hash_map cofaces(simplex, dim, n, modulus, dist, binomial_coeff); -// simplex_coboundary_enumerator cofaces_full(simplex, dim, n, modulus, dist_full, binomial_coeff); while (cofaces.has_next()) { -// assert(cofaces_full.has_next()); diameter_entry_t coface = cofaces.next(); -// , coface_full = cofaces_full.next(); -// value_t diameter = comp.diameter(get_index(coface)); -// assert(comp.diameter(get_index(coface)) == get_diameter(coface)); if (get_diameter(coface) <= threshold) { coface_entries.push_back(coface); @@ -1119,7 +1105,7 @@ int main(int argc, char** argv) { hash_map pivot_column_index; pivot_column_index.reserve(columns_to_reduce.size()); - compute_pairs(columns_to_reduce, pivot_column_index, dim, n, threshold, modulus, multiplicative_inverse, sparse_dist, //dist, + compute_pairs(columns_to_reduce, pivot_column_index, dim, n, threshold, modulus, multiplicative_inverse, sparse_dist, comp, comp_prev, binomial_coeff); if (dim < dim_max) { -- cgit v1.2.3 From 5adf0472e278034d155395cb8a95c017b30f0acb Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 29 Sep 2016 18:16:54 +0200 Subject: rearranged reduction columns code --- ripser.cpp | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index ec73f6d..4067018 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -685,7 +685,7 @@ void compute_pairs(std::vector& columns_to_reduce, hash_map reduction_matrix; + compressed_sparse_matrix reduction_coefficients; #else #ifdef USE_COEFFICIENTS std::vector reduction_coefficients; @@ -723,9 +723,9 @@ void compute_pairs(std::vector& columns_to_reduce, hash_map& columns_to_reduce, hash_map& columns_to_reduce, hash_map& columns_to_reduce, hash_map 0); #endif - reduction_matrix.push_back(e); + reduction_coefficients.push_back(e); } #else #ifdef USE_COEFFICIENTS -- cgit v1.2.3 From 43d28e6cea16062d0cefaf16979d07837cc788d6 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 29 Sep 2016 23:02:29 +0200 Subject: removed comparator --- ripser.cpp | 88 +++----------------------------------------------------------- 1 file changed, 3 insertions(+), 85 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 713c70b..f5384e9 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -213,42 +213,6 @@ template struct greater_diameter_or_smaller_index { } }; -template class rips_filtration_comparator { -public: - const DistanceMatrix& dist; - const index_t dim; - -private: - mutable std::vector vertices; - const binomial_coeff_table& binomial_coeff; - -public: - rips_filtration_comparator(const DistanceMatrix& _dist, const index_t _dim, - const binomial_coeff_table& _binomial_coeff) - : dist(_dist), dim(_dim), vertices(_dim + 1), binomial_coeff(_binomial_coeff){}; - - value_t diameter(const index_t index) const { - value_t diam = 0; - get_simplex_vertices(index, dim, dist.size(), binomial_coeff, vertices.begin()); - - 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])); } - return diam; - } - - bool operator()(const index_t a, const index_t b) const { - assert(a < binomial_coeff(dist.size(), dim + 1)); - assert(b < binomial_coeff(dist.size(), dim + 1)); - - return greater_diameter_or_smaller_index()(diameter_index_t(diameter(a), a), - diameter_index_t(diameter(b), b)); - } - - template bool operator()(const Entry& a, const Entry& b) const { - return operator()(get_index(a), get_index(b)); - } -}; - template class simplex_coboundary_enumerator { private: index_t idx_below, idx_above, v, k; @@ -584,51 +548,10 @@ template void push_entry(Heap& column, index_t i, coefficient_t column.push(std::make_pair(diameter, e)); } -template -void assemble_columns_to_reduce(std::vector& columns_to_reduce, - hash_map& pivot_column_index, const Comparator& comp, index_t dim, - index_t n, value_t threshold, const binomial_coeff_table& binomial_coeff) { - index_t num_simplices = binomial_coeff(n, dim + 2); - - 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 = comp.diameter(index); - 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"; -#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 assemble_columns_to_reduce(std::vector& simplices, std::vector& columns_to_reduce, hash_map& pivot_column_index, const sparse_distance_matrix& sparse_dist, - const Comparator& comp, index_t dim, index_t n, value_t threshold, const coefficient_t modulus, const binomial_coeff_table& binomial_coeff) { @@ -668,12 +591,10 @@ void assemble_columns_to_reduce(std::vector& simplices, #endif } -template +template void compute_pairs(std::vector& columns_to_reduce, hash_map& pivot_column_index, index_t dim, index_t n, value_t threshold, coefficient_t modulus, const std::vector& multiplicative_inverse, const DistanceMatrix& dist, - const ComparatorCofaces& comp, const Comparator& comp_prev, const binomial_coeff_table& binomial_coeff) { #ifdef PRINT_PERSISTENCE_PAIRS @@ -1099,17 +1020,14 @@ int main(int argc, char** argv) { } for (index_t dim = 1; dim <= dim_max; ++dim) { - rips_filtration_comparator comp(dist, dim + 1, binomial_coeff); - rips_filtration_comparator comp_prev(dist, dim, binomial_coeff); - hash_map pivot_column_index; pivot_column_index.reserve(columns_to_reduce.size()); compute_pairs(columns_to_reduce, pivot_column_index, dim, n, threshold, modulus, multiplicative_inverse, sparse_dist, - comp, comp_prev, binomial_coeff); + binomial_coeff); if (dim < dim_max) { - assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, sparse_dist, comp, dim, n, threshold, modulus, binomial_coeff); + assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, sparse_dist, dim, n, threshold, modulus, binomial_coeff); } } } -- cgit v1.2.3 From 77deb431e86bab1b4d3716909fc061ac52f692f1 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 29 Sep 2016 23:03:12 +0200 Subject: sparse graph support for dim 0 --- ripser.cpp | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index f5384e9..561f7bf 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -982,17 +982,22 @@ int main(int argc, char** argv) { std::vector columns_to_reduce; std::vector simplices , &edges = simplices; - rips_filtration_comparator comp(dist, 1, binomial_coeff); - for (index_t index = binomial_coeff(n, 2); index-- > 0;) { - value_t diameter = comp.diameter(index); - if (diameter <= threshold) edges.push_back(std::make_pair(diameter, index)); + + for (index_t i = 0; i < n; ++i) { + simplex_coboundary_enumerator cofaces(i, 0, n, modulus, sparse_dist, binomial_coeff); + + while (cofaces.has_next(false)) { + auto coface = cofaces.next(); + value_t diameter = get_diameter(coface); + if (diameter <= threshold) edges.push_back(std::make_pair(diameter, get_index(coface))); + } + + std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index()); } { union_find dset(n); - 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 -- cgit v1.2.3 From 7cd0e22d554207e420d4dcedbe7e753bac64c3c3 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Fri, 30 Sep 2016 18:43:48 +0200 Subject: pretty-print --- ripser.cpp | 90 +++++++++++++++++++++++++++++--------------------------------- 1 file changed, 42 insertions(+), 48 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 6faa9c4..3effc85 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -18,7 +18,6 @@ with this program. If not, see . */ - //#define ASSEMBLE_REDUCTION_MATRIX //#define USE_COEFFICIENTS @@ -167,7 +166,7 @@ typedef index_t entry_t; const index_t get_index(entry_t i) { return i; } index_t get_coefficient(entry_t i) { return 1; } entry_t make_entry(index_t _index, coefficient_t _value) { return entry_t(_index); } -void set_coefficient(index_t& e, const coefficient_t c) { } +void set_coefficient(index_t& e, const coefficient_t c) {} #endif @@ -177,7 +176,7 @@ template struct smaller_index { bool operator()(const Entry& a, const Entry& b) { return get_index(a) < get_index(b); } }; -class diameter_index_t: public std::pair { +class diameter_index_t : public std::pair { public: diameter_index_t() : std::pair() {} diameter_index_t(std::pair p) : std::pair(p) {} @@ -284,7 +283,8 @@ public: value_t coface_diameter = get_diameter(simplex); for (index_t w : vertices) coface_diameter = std::max(coface_diameter, dist(v, w)); coefficient_t coface_coefficient = (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; - return diameter_entry_t(coface_diameter, idx_above + binomial_coeff(v--, k + 1) + idx_below, coface_coefficient); + return diameter_entry_t(coface_diameter, idx_above + binomial_coeff(v--, k + 1) + idx_below, + coface_coefficient); } }; @@ -340,9 +340,7 @@ template struct second_argument_equal_to { bool operator()(const T& lhs, const T& rhs) const { return lhs.second == rhs.second; } }; - -template<> -class simplex_coboundary_enumerator { +template <> class simplex_coboundary_enumerator { private: index_t idx_below, idx_above, v, k, max_vertex_below; const binomial_coeff_table& binomial_coeff; @@ -359,8 +357,8 @@ private: public: simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, index_t _n, - const coefficient_t _modulus, const sparse_distance_matrix& _sparse_dist, - const binomial_coeff_table& _binomial_coeff) + const coefficient_t _modulus, const sparse_distance_matrix& _sparse_dist, + const binomial_coeff_table& _binomial_coeff) : simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(_n - 1), k(_dim + 1), max_vertex_below(_n - 1), modulus(_modulus), sparse_dist(_sparse_dist), binomial_coeff(_binomial_coeff) { get_simplex_vertices(idx_below, _dim, _n, _binomial_coeff, std::back_inserter(vertices)); @@ -384,7 +382,8 @@ public: else x = std::min(x, y, greater_diameter_or_smaller_index()); } - return all_cofaces || !(k > 0 && get_next_vertex(max_vertex_below, idx_below, k, binomial_coeff) > get_index(x)); + return all_cofaces || + !(k > 0 && get_next_vertex(max_vertex_below, idx_below, k, binomial_coeff) > get_index(x)); continue_outer:; } return false; @@ -403,9 +402,8 @@ public: 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, idx_above + binomial_coeff(get_index(x), k + 1) + idx_below, + coface_coefficient); } }; @@ -595,9 +593,10 @@ void assemble_columns_to_reduce(std::vector& columns_to_reduce value_t diameter = comp.diameter(index); 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 } } @@ -614,16 +613,13 @@ void assemble_columns_to_reduce(std::vector& columns_to_reduce #endif } - - template void assemble_columns_to_reduce(std::vector& simplices, - std::vector& columns_to_reduce, - hash_map& pivot_column_index, - const sparse_distance_matrix& sparse_dist, - const Comparator& comp, - index_t dim, index_t n, value_t threshold, const coefficient_t modulus, - const binomial_coeff_table& binomial_coeff) { + std::vector& columns_to_reduce, + hash_map& pivot_column_index, + const sparse_distance_matrix& sparse_dist, const Comparator& comp, index_t dim, + index_t n, value_t threshold, const coefficient_t modulus, + const binomial_coeff_table& binomial_coeff) { #ifdef INDICATE_PROGRESS std::cout << "\033[K" @@ -631,24 +627,25 @@ void assemble_columns_to_reduce(std::vector& simplices, #endif columns_to_reduce.clear(); - + std::vector next_simplices; for (diameter_index_t simplex : simplices) { - simplex_coboundary_enumerator cofaces(simplex, dim, n, modulus, sparse_dist, binomial_coeff); + simplex_coboundary_enumerator cofaces(simplex, dim, n, modulus, sparse_dist, + binomial_coeff); 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"; @@ -661,12 +658,10 @@ void assemble_columns_to_reduce(std::vector& simplices, #endif } -template +template void compute_pairs(std::vector& columns_to_reduce, hash_map& pivot_column_index, index_t dim, index_t n, value_t threshold, coefficient_t modulus, - const std::vector& multiplicative_inverse, - const DistanceMatrix& dist, + const std::vector& multiplicative_inverse, const DistanceMatrix& dist, const ComparatorCofaces& comp, const Comparator& comp_prev, const binomial_coeff_table& binomial_coeff) { @@ -736,8 +731,7 @@ void compute_pairs(std::vector& columns_to_reduce, hash_map columns_to_reduce; - std::vector simplices , &edges = simplices; - rips_filtration_comparator comp(dist, 1, binomial_coeff); - for (index_t index = binomial_coeff(n, 2); index-- > 0;) { - value_t diameter = comp.diameter(index); - if (diameter <= threshold) edges.push_back(std::make_pair(diameter, index)); - } - + std::vector simplices, &edges = simplices; + rips_filtration_comparator comp(dist, 1, binomial_coeff); + for (index_t index = binomial_coeff(n, 2); index-- > 0;) { + value_t diameter = comp.diameter(index); + if (diameter <= threshold) edges.push_back(std::make_pair(diameter, index)); + } + { union_find dset(n); @@ -1071,8 +1065,7 @@ int main(int argc, char** argv) { if (u != v) { #ifdef PRINT_PERSISTENCE_PAIRS - if (get_diameter(e) > 0) - std::cout << " [0," << get_diameter(e) << ")" << std::endl; + if (get_diameter(e) > 0) std::cout << " [0," << get_diameter(e) << ")" << std::endl; #endif dset.link(u, v); } else @@ -1093,11 +1086,12 @@ int main(int argc, char** argv) { hash_map pivot_column_index; pivot_column_index.reserve(columns_to_reduce.size()); - compute_pairs(columns_to_reduce, pivot_column_index, dim, n, threshold, modulus, multiplicative_inverse, sparse_dist, - comp, comp_prev, binomial_coeff); + compute_pairs(columns_to_reduce, pivot_column_index, dim, n, threshold, modulus, multiplicative_inverse, + sparse_dist, comp, comp_prev, binomial_coeff); if (dim < dim_max) { - assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, sparse_dist, comp, dim, n, threshold, modulus, binomial_coeff); + assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, sparse_dist, comp, dim, n, + threshold, modulus, binomial_coeff); } } } -- cgit v1.2.3 From 7f2d16e1054529bed4d04c00a6fd7a7c33fdd458 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Fri, 30 Sep 2016 22:03:08 +0200 Subject: fixed build problems --- ripser.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index b183b3d..36a477a 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -180,6 +180,7 @@ class diameter_index_t : public std::pair { public: diameter_index_t() : std::pair() {} diameter_index_t(std::pair p) : std::pair(p) {} + diameter_index_t(index_t i) : std::pair(0, i) {} }; value_t get_diameter(diameter_index_t i) { return i.first; } index_t get_index(diameter_index_t i) { return i.second; } @@ -187,8 +188,6 @@ index_t get_index(diameter_index_t i) { return i.second; } class diameter_entry_t : public std::pair { public: diameter_entry_t(std::pair p) : std::pair(p) {} - diameter_entry_t(entry_t e) : std::pair(0, e) {} - diameter_entry_t() : diameter_entry_t(0) {} diameter_entry_t(value_t _diameter, index_t _index, coefficient_t _coefficient) : std::pair(_diameter, make_entry(_index, _coefficient)) {} diameter_entry_t(diameter_index_t _diameter_index, coefficient_t _coefficient) @@ -968,7 +967,7 @@ int main(int argc, char** argv) { std::vector simplices, &edges = simplices; for (index_t i = 0; i < n; ++i) { - simplex_coboundary_enumerator cofaces(i, 0, n, modulus, sparse_dist, binomial_coeff); + simplex_coboundary_enumerator cofaces(diameter_index_t(i), 0, n, modulus, sparse_dist, binomial_coeff); while (cofaces.has_next(false)) { auto coface = cofaces.next(); -- cgit v1.2.3 From 0e3bd6152519cce4a55c104aa00478aeaf23dfb9 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 9 Nov 2016 23:58:52 +0100 Subject: initialize edges from lower distance matrix --- ripser.cpp | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 36a477a..3ebf2e3 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -965,19 +965,19 @@ int main(int argc, char** argv) { std::vector columns_to_reduce; std::vector simplices, &edges = simplices; - - for (index_t i = 0; i < n; ++i) { - simplex_coboundary_enumerator cofaces(diameter_index_t(i), 0, n, modulus, sparse_dist, binomial_coeff); - - while (cofaces.has_next(false)) { - auto coface = cofaces.next(); - value_t diameter = get_diameter(coface); - if (diameter <= threshold) edges.push_back(std::make_pair(diameter, get_index(coface))); - } - - std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index()); + + std::chrono::time_point start, end; + long elapsed_time; + + start = std::chrono::system_clock::now(); + + for (index_t e = 0; e < dist.distances.size(); ++e) { + value_t diameter = dist.distances[e]; + if (diameter <= threshold) edges.push_back(std::make_pair(diameter, e)); } - + + std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index()); + { union_find dset(n); @@ -1019,4 +1019,9 @@ int main(int argc, char** argv) { threshold, modulus, binomial_coeff); } } + + end = std::chrono::system_clock::now(); + elapsed_time = std::chrono::duration_cast(end - start).count(); + std::cout << "Computed Rips persistence in " << elapsed_time << " ms.\n"; + } -- cgit v1.2.3 From 9213750960682a4bd3d4291467818395f124355c Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 10 Nov 2016 00:25:10 +0100 Subject: more timing --- ripser.cpp | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 3ebf2e3..d984904 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -946,8 +946,14 @@ int main(int argc, char** argv) { exit(-1); } + std::chrono::time_point start; + + start = std::chrono::system_clock::now(); + compressed_lower_distance_matrix dist = read_file(filename ? file_stream : std::cin, format); + std::cout << "Reading file in " << std::chrono::duration_cast(std::chrono::system_clock::now() - start).count() / 1000. << " s.\n"; + index_t n = dist.size(); std::cout << "distance matrix with " << n << " points" << std::endl; @@ -955,7 +961,12 @@ int main(int argc, char** argv) { auto value_range = std::minmax_element(dist.distances.begin(), dist.distances.end()); std::cout << "value range: [" << *value_range.first << "," << *value_range.second << "]" << std::endl; + start = std::chrono::system_clock::now(); + sparse_distance_matrix sparse_dist(dist, threshold); + + std::cout << "Building sparse distance matrix in " << std::chrono::duration_cast(std::chrono::system_clock::now() - start).count() / 1000. << " s.\n"; + dim_max = std::min(dim_max, n - 2); @@ -966,8 +977,6 @@ int main(int argc, char** argv) { std::vector simplices, &edges = simplices; - std::chrono::time_point start, end; - long elapsed_time; start = std::chrono::system_clock::now(); @@ -1020,8 +1029,7 @@ int main(int argc, char** argv) { } } - end = std::chrono::system_clock::now(); - elapsed_time = std::chrono::duration_cast(end - start).count(); - std::cout << "Computed Rips persistence in " << elapsed_time << " ms.\n"; + std::cout << "Computing Rips persistence in " << std::chrono::duration_cast(std::chrono::system_clock::now() - start).count() / 1000. << " s.\n"; + } -- cgit v1.2.3 From c86c7fc1620363b4a78699c88644d5d8c6051379 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 30 Nov 2016 17:46:00 -0500 Subject: cleanup --- ripser.cpp | 30 ++++++++++++------------------ 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 270f8b9..a66e8e6 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -137,6 +137,8 @@ index_t get_index(diameter_index_t i) { return i.second; } class diameter_entry_t : public std::pair { public: diameter_entry_t(std::pair p) : std::pair(p) {} + diameter_entry_t(entry_t e) : std::pair(0, e) {} + diameter_entry_t() : diameter_entry_t(0) {} diameter_entry_t(value_t _diameter, index_t _index, coefficient_t _coefficient) : std::pair(_diameter, make_entry(_index, _coefficient)) {} diameter_entry_t(diameter_index_t _diameter_index, coefficient_t _coefficient) @@ -203,14 +205,6 @@ public: size_t size() const { return neighbors.size(); } }; -template struct second_argument_greater { - bool operator()(const T& lhs, const T& rhs) const { return lhs.second > rhs.second; } -}; - -template struct second_argument_equal_to { - bool operator()(const T& lhs, const T& rhs) const { return lhs.second == rhs.second; } -}; - template <> void compressed_distance_matrix::init_rows() { value_t* pointer = &distances[0]; for (index_t i = 1; i < size(); ++i) { @@ -494,16 +488,16 @@ public: return out; } -// value_t compute_diameter(const index_t index, index_t dim) const { -// value_t diam = -std::numeric_limits::infinity(); -// -// vertices.clear(); -// get_simplex_vertices(index, dim, sparse_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, sparse_dist(vertices[i], vertices[j])); } -// return diam; -// } + value_t compute_diameter(const index_t index, index_t dim) const { + value_t diam = -std::numeric_limits::infinity(); + + vertices.clear(); + get_simplex_vertices(index, dim, sparse_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])); } + return diam; + } void assemble_columns_to_reduce(std::vector& simplices, std::vector& columns_to_reduce, -- cgit v1.2.3 From 117088b16d9abe1ca18a909568bd71e377782ff4 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 30 Nov 2016 18:22:01 -0500 Subject: speedup coboundary enumeration --- ripser.cpp | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index a66e8e6..aa2d438 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -381,7 +381,10 @@ class ripser { coefficient_t modulus; compressed_lower_distance_matrix dist; sparse_distance_matrix sparse_dist; + mutable std::vector vertices; + mutable std::vector::const_reverse_iterator> ii; + mutable std::vector::const_reverse_iterator> ee; public: ripser(compressed_lower_distance_matrix&& _dist,sparse_distance_matrix&& _sparse_dist, index_t _dim_max, value_t _threshold, coefficient_t _modulus) @@ -400,9 +403,11 @@ public: const sparse_distance_matrix& sparse_dist; const binomial_coeff_table& binomial_coeff; - std::vector vertices; + std::vector& vertices; - std::vector::const_reverse_iterator> ii, ee; + std::vector::const_reverse_iterator>& ii; + std::vector::const_reverse_iterator>& ee; + diameter_index_t x; @@ -410,7 +415,12 @@ public: public: simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, const ripser& _parent) : parent(_parent), simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), k(_dim + 1), - max_vertex_below(parent.n - 1), modulus(parent.modulus), sparse_dist(parent.sparse_dist), binomial_coeff(parent.binomial_coeff) { + max_vertex_below(parent.n - 1), modulus(parent.modulus), sparse_dist(parent.sparse_dist), binomial_coeff(parent.binomial_coeff), vertices(parent.vertices), ii(parent.ii), ee(parent.ee) { + + ii.clear(); + ee.clear(); + vertices.clear(); + parent.get_simplex_vertices(idx_below, _dim, parent.n, std::back_inserter(vertices)); for (auto v : vertices) { -- cgit v1.2.3 From 63c3ae66f90bbf4ef295433b02d8dd10072bcf6d Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 30 Nov 2016 19:23:17 -0500 Subject: use only sparse distance matrix --- ripser.cpp | 156 +++++++++++++++++++++++++++---------------------------------- 1 file changed, 68 insertions(+), 88 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index aa2d438..01e587b 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,7 +88,6 @@ 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)); @@ -379,22 +379,20 @@ class ripser { const binomial_coeff_table binomial_coeff; std::vector multiplicative_inverse; coefficient_t modulus; - compressed_lower_distance_matrix dist; sparse_distance_matrix sparse_dist; - + mutable std::vector vertices; mutable std::vector::const_reverse_iterator> ii; mutable std::vector::const_reverse_iterator> ee; public: - ripser(compressed_lower_distance_matrix&& _dist,sparse_distance_matrix&& _sparse_dist, index_t _dim_max, value_t _threshold, coefficient_t _modulus) - : dist(_dist), sparse_dist(_sparse_dist), n(_sparse_dist.size()), dim_max(std::min(_dim_max, index_t(_sparse_dist.size() - 2))), threshold(_threshold), - modulus(_modulus), binomial_coeff(n, dim_max + 2), - multiplicative_inverse(multiplicative_inverse_vector(_modulus)) {} - + ripser(sparse_distance_matrix&& _sparse_dist, index_t _dim_max, value_t _threshold, coefficient_t _modulus) + : sparse_dist(_sparse_dist), n(_sparse_dist.size()), + dim_max(std::min(_dim_max, index_t(_sparse_dist.size() - 2))), threshold(_threshold), modulus(_modulus), + binomial_coeff(n, dim_max + 2), multiplicative_inverse(multiplicative_inverse_vector(_modulus)) {} + class simplex_coboundary_enumerator { private: - const ripser& parent; index_t idx_below, idx_above, v, k, max_vertex_below; @@ -404,31 +402,28 @@ public: const binomial_coeff_table& binomial_coeff; std::vector& vertices; - std::vector::const_reverse_iterator>& ii; std::vector::const_reverse_iterator>& ee; - - diameter_index_t x; - - + public: simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, const ripser& _parent) - : parent(_parent), simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), k(_dim + 1), - max_vertex_below(parent.n - 1), modulus(parent.modulus), sparse_dist(parent.sparse_dist), binomial_coeff(parent.binomial_coeff), vertices(parent.vertices), ii(parent.ii), ee(parent.ee) { - + : parent(_parent), simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), + k(_dim + 1), max_vertex_below(parent.n - 1), modulus(parent.modulus), sparse_dist(parent.sparse_dist), + binomial_coeff(parent.binomial_coeff), vertices(parent.vertices), ii(parent.ii), ee(parent.ee) { + ii.clear(); ee.clear(); vertices.clear(); - + parent.get_simplex_vertices(idx_below, _dim, parent.n, std::back_inserter(vertices)); - + for (auto v : vertices) { ii.push_back(sparse_dist.neighbors[v].rbegin()); ee.push_back(sparse_dist.neighbors[v].rend()); } } - + bool has_next(bool all_cofaces = true) { for (auto &it0 = ii[0], &end0 = ee[0]; it0 != end0; ++it0) { x = *it0; @@ -440,34 +435,32 @@ public: if (get_index(y) != get_index(x)) goto continue_outer; else - x = std::min(x, y, greater_diameter_or_smaller_index()); + x = std::max(x, y); } - return all_cofaces || - !(k > 0 && parent.get_next_vertex(max_vertex_below, idx_below, k) > get_index(x)); + return all_cofaces || !(k > 0 && parent.get_next_vertex(max_vertex_below, idx_below, k) > get_index(x)); continue_outer:; } return false; } - + diameter_entry_t next() { ++ii[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)); - + 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); + coface_coefficient); } }; - 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; @@ -485,10 +478,11 @@ public: assert(binomial_coeff(v, k) <= idx && binomial_coeff(v + 1, k) > idx); return v; } - + + index_t get_edge_index(const index_t i, const index_t j) const { return binomial_coeff(i, 2) + j; } + template - OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v, - OutputIterator out) const { + OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v, OutputIterator out) const { --v; for (index_t k = dim + 1; k > 0; --k) { get_next_vertex(v, idx, k); @@ -497,57 +491,46 @@ public: } return out; } - - value_t compute_diameter(const index_t index, index_t dim) const { - value_t diam = -std::numeric_limits::infinity(); - vertices.clear(); - get_simplex_vertices(index, dim, sparse_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])); } - return diam; - } - -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, + hash_map& pivot_column_index, index_t dim) { #ifdef INDICATE_PROGRESS - std::cout << "\033[K" - << "assembling columns" << std::flush << "\r"; + std::cout << "\033[K" + << "assembling columns" << std::flush << "\r"; #endif - columns_to_reduce.clear(); + columns_to_reduce.clear(); - std::vector next_simplices; + std::vector next_simplices; - for (diameter_index_t simplex : simplices) { - simplex_coboundary_enumerator cofaces(simplex, dim, *this); + for (diameter_index_t simplex : simplices) { + simplex_coboundary_enumerator cofaces(simplex, dim, *this); - while (cofaces.has_next(false)) { - auto coface = cofaces.next(); + while (cofaces.has_next(false)) { + auto coface = cofaces.next(); - next_simplices.push_back(std::make_pair(get_diameter(coface), get_index(coface))); + 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))); + 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); + simplices.swap(next_simplices); #ifdef INDICATE_PROGRESS - std::cout << "\033[K" - << "sorting " << columns_to_reduce.size() << " columns" << std::flush << "\r"; + 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()); + 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) { @@ -723,29 +706,30 @@ void assemble_columns_to_reduce(std::vector& simplices, void compute_barcodes() { std::vector columns_to_reduce; - + std::vector simplices, &edges = simplices; - - for (index_t e = 0; e < dist.distances.size(); ++e) { - value_t diameter = dist.distances[e]; - if (diameter <= threshold) edges.push_back(std::make_pair(diameter, e)); - } - + + for (index_t i = 0; i < n; ++i) + for (auto n : sparse_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))); + } + std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index()); - + { union_find dset(n); - + #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; @@ -755,24 +739,21 @@ void assemble_columns_to_reduce(std::vector& simplices, 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(simplices, columns_to_reduce, pivot_column_index, dim); - } + + if (dim < dim_max) { assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, dim); } } - } }; @@ -995,9 +976,8 @@ int main(int argc, char** argv) { auto value_range = std::minmax_element(dist.distances.begin(), dist.distances.end()); std::cout << "value range: [" << *value_range.first << "," << *value_range.second << "]" << std::endl; - + sparse_distance_matrix sparse_dist(dist, threshold); - - ripser(std::move(dist), std::move(sparse_dist), dim_max, threshold, modulus).compute_barcodes(); + ripser(std::move(sparse_dist), dim_max, threshold, modulus).compute_barcodes(); } -- cgit v1.2.3 From f1af637b5d3552472e63a6715ea073540871bf65 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 28 Dec 2016 14:32:21 +0100 Subject: code cleanup --- ripser.cpp | 126 ++++++++++++++++++++++++++++++++----------------------------- 1 file changed, 66 insertions(+), 60 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 01e587b..ec06d05 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -24,6 +24,8 @@ with this program. If not, see . //#define INDICATE_PROGRESS #define PRINT_PERSISTENCE_PAIRS +#define SPARSE_DISTANCE_MATRIX + //#define USE_GOOGLE_HASHMAP #include @@ -52,7 +54,6 @@ 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++) { @@ -129,7 +130,6 @@ class diameter_index_t : public std::pair { public: diameter_index_t() : std::pair() {} diameter_index_t(std::pair p) : std::pair(p) {} - diameter_index_t(index_t i) : std::pair(0, i) {} }; value_t get_diameter(diameter_index_t i) { return i.first; } index_t get_index(diameter_index_t i) { return i.second; } @@ -379,51 +379,88 @@ class ripser { const binomial_coeff_table binomial_coeff; std::vector multiplicative_inverse; coefficient_t modulus; +#ifdef SPARSE_DISTANCE_MATRIX sparse_distance_matrix sparse_dist; - +#else + compressed_lower_distance_matrix dist; +#endif mutable std::vector vertices; +#ifdef SPARSE_DISTANCE_MATRIX mutable std::vector::const_reverse_iterator> ii; mutable std::vector::const_reverse_iterator> ee; - +#endif + public: ripser(sparse_distance_matrix&& _sparse_dist, index_t _dim_max, value_t _threshold, coefficient_t _modulus) : sparse_dist(_sparse_dist), n(_sparse_dist.size()), dim_max(std::min(_dim_max, index_t(_sparse_dist.size() - 2))), 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; + while (count > 0) { + index_t i = v; + index_t step = count >> 1; + i -= step; + if (binomial_coeff(i, k) > idx) { + v = --i; + count -= step + 1; + } else + count = step; + } + } + assert(binomial_coeff(v, k) <= idx && binomial_coeff(v + 1, k) > idx); + return v; + } + + index_t get_edge_index(const index_t i, const index_t j) const { return binomial_coeff(i, 2) + j; } + + template + OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v, + OutputIterator out) const { + --v; + for (index_t k = dim + 1; k > 0; --k) { + get_next_vertex(v, idx, k); + *out++ = v; + idx -= binomial_coeff(v, k); + } + return out; + } + class simplex_coboundary_enumerator { private: const ripser& parent; - + index_t idx_below, idx_above, v, k, max_vertex_below; const diameter_entry_t simplex; const coefficient_t modulus; const sparse_distance_matrix& sparse_dist; const binomial_coeff_table& binomial_coeff; - + std::vector& vertices; std::vector::const_reverse_iterator>& ii; std::vector::const_reverse_iterator>& ee; diameter_index_t x; - + public: simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, const ripser& _parent) - : parent(_parent), simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), - k(_dim + 1), max_vertex_below(parent.n - 1), modulus(parent.modulus), sparse_dist(parent.sparse_dist), - binomial_coeff(parent.binomial_coeff), vertices(parent.vertices), ii(parent.ii), ee(parent.ee) { - + : parent(_parent), simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), + k(_dim + 1), max_vertex_below(parent.n - 1), modulus(parent.modulus), sparse_dist(parent.sparse_dist), + binomial_coeff(parent.binomial_coeff), vertices(parent.vertices), ii(parent.ii), ee(parent.ee) { + ii.clear(); ee.clear(); vertices.clear(); - + parent.get_simplex_vertices(idx_below, _dim, parent.n, std::back_inserter(vertices)); - + for (auto v : vertices) { ii.push_back(sparse_dist.neighbors[v].rbegin()); ee.push_back(sparse_dist.neighbors[v].rend()); } } - + bool has_next(bool all_cofaces = true) { for (auto &it0 = ii[0], &end0 = ee[0]; it0 != end0; ++it0) { x = *it0; @@ -442,55 +479,25 @@ public: } return false; } - + diameter_entry_t next() { ++ii[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)); - + 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); + coface_coefficient); } }; - 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; - while (count > 0) { - index_t i = v; - index_t step = count >> 1; - i -= step; - if (binomial_coeff(i, k) > idx) { - v = --i; - count -= step + 1; - } else - count = step; - } - } - assert(binomial_coeff(v, k) <= idx && binomial_coeff(v + 1, k) > idx); - return v; - } - - index_t get_edge_index(const index_t i, const index_t j) const { return binomial_coeff(i, 2) + j; } - - template - OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v, OutputIterator out) const { - --v; - for (index_t k = dim + 1; k > 0; --k) { - get_next_vertex(v, idx, k); - *out++ = v; - idx -= binomial_coeff(v, k); - } - return out; - } void assemble_columns_to_reduce(std::vector& simplices, std::vector& columns_to_reduce, @@ -707,18 +714,17 @@ public: std::vector columns_to_reduce; - std::vector simplices, &edges = simplices; - - for (index_t i = 0; i < n; ++i) - for (auto n : sparse_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))); - } - - std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index()); - + std::vector simplices; + { union_find dset(n); + std::vector &edges = simplices; + for (index_t i = 0; i < n; ++i) + for (auto n : sparse_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))); + } + std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index()); #ifdef PRINT_PERSISTENCE_PAIRS std::cout << "persistence intervals in dim 0:" << std::endl; -- cgit v1.2.3 From 31d0819c553fabbb25192797b250967216d699d3 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 29 Dec 2016 16:12:38 +0100 Subject: more code cleanup --- ripser.cpp | 317 +++++++++++++++++++++++++++++++++++++++++++------------------ 1 file changed, 224 insertions(+), 93 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index ec06d05..d0674f9 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -221,14 +221,12 @@ template <> void compressed_distance_matrix::init_rows() { } } -template <> value_t compressed_distance_matrix::operator()(index_t i, index_t j) const { - if (i > j) std::swap(i, j); - return i == j ? 0 : 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]; } -template <> value_t compressed_distance_matrix::operator()(index_t i, index_t j) const { - if (i > j) std::swap(i, j); - return i == j ? 0 : rows[j][i]; +template <> value_t compressed_distance_matrix::operator()(const index_t i, const index_t j) const { + return i == j ? 0 : i > j ? rows[i][j] : rows[j][i]; } typedef compressed_distance_matrix compressed_lower_distance_matrix; @@ -373,28 +371,22 @@ template void push_entry(Heap& column, index_t i, coefficient_t column.push(std::make_pair(diameter, e)); } -class ripser { - index_t dim_max, n; +template class ripser { + DistanceMatrix dist; + index_t n, dim_max; value_t threshold; const binomial_coeff_table binomial_coeff; std::vector multiplicative_inverse; coefficient_t modulus; -#ifdef SPARSE_DISTANCE_MATRIX - sparse_distance_matrix sparse_dist; -#else - compressed_lower_distance_matrix dist; -#endif mutable std::vector vertices; -#ifdef SPARSE_DISTANCE_MATRIX - mutable std::vector::const_reverse_iterator> ii; - mutable std::vector::const_reverse_iterator> ee; -#endif - + mutable std::vector::const_reverse_iterator> neighbor_it; + mutable std::vector::const_reverse_iterator> neighbor_end; + public: - ripser(sparse_distance_matrix&& _sparse_dist, index_t _dim_max, value_t _threshold, coefficient_t _modulus) - : sparse_dist(_sparse_dist), n(_sparse_dist.size()), - dim_max(std::min(_dim_max, index_t(_sparse_dist.size() - 2))), threshold(_threshold), modulus(_modulus), - binomial_coeff(n, dim_max + 2), multiplicative_inverse(multiplicative_inverse_vector(_modulus)) {} + ripser(DistanceMatrix&& _dist, index_t _dim_max, value_t _threshold, coefficient_t _modulus) + : dist(std::move(_dist)), n(dist.size()), dim_max(std::min(_dim_max, index_t(dist.size() - 2))), + 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) { @@ -417,8 +409,7 @@ public: index_t get_edge_index(const index_t i, const index_t j) const { return binomial_coeff(i, 2) + j; } template - OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v, - OutputIterator out) const { + OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v, OutputIterator out) const { --v; for (index_t k = dim + 1; k > 0; --k) { get_next_vertex(v, idx, k); @@ -427,45 +418,92 @@ public: } return out; } - + + value_t compute_diameter(const index_t index, index_t dim) const { + value_t diam = -std::numeric_limits::infinity(); + + vertices.clear(); + 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])); } + 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) + : simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), k(_dim + 1), + modulus(parent.modulus), binomial_coeff(parent.binomial_coeff), dist(parent.dist), vertices(_dim + 1) { + 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); + } + }; + + class simplex_sparse_coboundary_enumerator { private: const ripser& parent; - + index_t idx_below, idx_above, v, k, max_vertex_below; const diameter_entry_t simplex; const coefficient_t modulus; - const sparse_distance_matrix& sparse_dist; + const DistanceMatrix& dist; const binomial_coeff_table& binomial_coeff; - + std::vector& vertices; - std::vector::const_reverse_iterator>& ii; - std::vector::const_reverse_iterator>& ee; + std::vector::const_reverse_iterator>& neighbor_it; + std::vector::const_reverse_iterator>& neighbor_end; diameter_index_t x; - + public: - simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, const ripser& _parent) - : parent(_parent), simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), - k(_dim + 1), max_vertex_below(parent.n - 1), modulus(parent.modulus), sparse_dist(parent.sparse_dist), - binomial_coeff(parent.binomial_coeff), vertices(parent.vertices), ii(parent.ii), ee(parent.ee) { - - ii.clear(); - ee.clear(); + simplex_sparse_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, const ripser& _parent) + : parent(_parent), simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), + k(_dim + 1), max_vertex_below(parent.n - 1), modulus(parent.modulus), dist(parent.dist), + binomial_coeff(parent.binomial_coeff), vertices(parent.vertices), neighbor_it(parent.neighbor_it), neighbor_end(parent.neighbor_end) { + + neighbor_it.clear(); + neighbor_end.clear(); vertices.clear(); - + parent.get_simplex_vertices(idx_below, _dim, parent.n, std::back_inserter(vertices)); - + for (auto v : vertices) { - ii.push_back(sparse_dist.neighbors[v].rbegin()); - ee.push_back(sparse_dist.neighbors[v].rend()); + neighbor_it.push_back(dist.neighbors[v].rbegin()); + neighbor_end.push_back(dist.neighbors[v].rend()); } } - + bool has_next(bool all_cofaces = true) { - for (auto &it0 = ii[0], &end0 = ee[0]; it0 != end0; ++it0) { + for (auto &it0 = neighbor_it[0], &end0 = neighbor_end[0]; it0 != end0; ++it0) { x = *it0; - for (size_t idx = 1; idx < ii.size(); ++idx) { - auto &it = ii[idx], end = ee[idx]; + 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)) if (++it == end) return false; auto y = *it; @@ -479,41 +517,77 @@ public: } return false; } - + diameter_entry_t next() { - ++ii[0]; - + ++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)); - + 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); + coface_coefficient); } }; - - void assemble_columns_to_reduce(std::vector& simplices, - std::vector& columns_to_reduce, + 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(); + +#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) % 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"; +#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 assemble_sparse_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(simplex, dim, *this); + simplex_sparse_coboundary_enumerator cofaces(simplex, dim, *this); while (cofaces.has_next(false)) { auto coface = cofaces.next(); @@ -539,6 +613,7 @@ public: #endif } + template void compute_pairs(std::vector& columns_to_reduce, hash_map& pivot_column_index, index_t dim) { @@ -624,7 +699,7 @@ public: #endif coface_entries.clear(); - simplex_coboundary_enumerator cofaces(simplex, dim, *this); + BoundaryEnumerator cofaces(simplex, dim, *this); while (cofaces.has_next()) { diameter_entry_t coface = cofaces.next(); if (get_diameter(coface) <= threshold) { @@ -710,58 +785,111 @@ public: #endif } - void compute_barcodes() { + void compute_barcodes(); +}; - std::vector columns_to_reduce; +template <> void ripser::compute_barcodes() { - std::vector simplices; - - { - union_find dset(n); - std::vector &edges = simplices; - for (index_t i = 0; i < n; ++i) - for (auto n : sparse_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))); - } - std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index()); + 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; + 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]); + 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) { + if (u != v) { #ifdef PRINT_PERSISTENCE_PAIRS - if (get_diameter(e) > 0) std::cout << " [0," << get_diameter(e) << ")" << std::endl; + 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); + 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); } + } +} + +template <> void ripser::compute_barcodes() { + + std::vector columns_to_reduce; + + std::vector simplices; + + { + union_find dset(n); + std::vector& edges = simplices; + 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))); } - std::reverse(columns_to_reduce.begin(), columns_to_reduce.end()); + 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 - for (index_t i = 0; i < n; ++i) - if (dset.find(i) == i) std::cout << " [0, )" << std::endl << std::flush; + 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()); + 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); + 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); } + if (dim < dim_max) { + assemble_sparse_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, dim + 1); } } -}; +} enum file_format { LOWER_DISTANCE_MATRIX, UPPER_DISTANCE_MATRIX, DISTANCE_MATRIX, POINT_CLOUD, DIPHA, RIPSER }; @@ -983,7 +1111,10 @@ int main(int argc, char** argv) { auto value_range = std::minmax_element(dist.distances.begin(), dist.distances.end()); std::cout << "value range: [" << *value_range.first << "," << *value_range.second << "]" << std::endl; - sparse_distance_matrix sparse_dist(dist, threshold); - ripser(std::move(sparse_dist), dim_max, threshold, modulus).compute_barcodes(); +#ifdef SPARSE_DISTANCE_MATRIX + ripser(std::move(dist), dim_max, threshold, modulus).compute_barcodes(); +#else + ripser(sparse_distance_matrix(dist, threshold), dim_max, threshold, modulus).compute_barcodes(); +#endif } -- cgit v1.2.3 From b945df6609d9198c34a982e70dee0773ee630897 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 29 Dec 2016 19:26:55 +0100 Subject: reordered initializers --- ripser.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 49ff6a8..3054441 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -508,8 +508,8 @@ public: public: simplex_sparse_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, const ripser& _parent) - : parent(_parent), simplex(_simplex), idx_below(get_index(_simplex)), idx_above(0), - v(parent.n - 1), k(_dim + 1), max_vertex_below(parent.n - 1), modulus(parent.modulus), + : 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) { -- cgit v1.2.3 From 2989ecd1047fcae2b8fead4b8fca1a4608342576 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 29 Dec 2016 23:07:00 +0100 Subject: use sparse distance matrix when threshold is given --- ripser.cpp | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 3054441..4ac5033 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -24,8 +24,6 @@ with this program. If not, see . //#define INDICATE_PROGRESS #define PRINT_PERSISTENCE_PAIRS -#define SPARSE_DISTANCE_MATRIX - //#define USE_GOOGLE_HASHMAP #include @@ -1063,14 +1061,14 @@ int main(int argc, char** argv) { std::cout << "value range: [" << *value_range.first << "," << *value_range.second << "]" << std::endl; -#ifdef SPARSE_DISTANCE_MATRIX - ripser(std::move(dist), dim_max, threshold, modulus) - .compute_barcodes(); -#else - ripser(sparse_distance_matrix(dist, threshold), dim_max, threshold, - modulus) - .compute_barcodes(); -#endif + if (threshold == std::numeric_limits::max()) + ripser(std::move(dist), dim_max, threshold, modulus) + .compute_barcodes(); + else + ripser(sparse_distance_matrix(dist, threshold), dim_max, threshold, + modulus) + .compute_barcodes(); + } template <> void ripser::compute_barcodes() { -- cgit v1.2.3 From 31a18897fb52e4e98faed5562e6b97da1e024b4c Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sat, 7 Jan 2017 00:31:25 +0100 Subject: move --- ripser.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 701ba35..fbbd345 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -1071,8 +1071,8 @@ int main(int argc, char** argv) { ripser(std::move(dist), dim_max, threshold, modulus) .compute_barcodes(); else - ripser(sparse_distance_matrix(dist, threshold), dim_max, threshold, - modulus) + ripser(sparse_distance_matrix(std::move(dist), threshold), dim_max, + threshold, modulus) .compute_barcodes(); } -- cgit v1.2.3 From 02e1d32f8c19a92b91dab3fc2520a6f4a98640f1 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 24 Aug 2017 12:36:23 +0200 Subject: threshold for sparse distance matrix --- ripser.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ripser.cpp b/ripser.cpp index 2f5a01a..74f9224 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -1079,7 +1079,7 @@ int main(int argc, char** argv) { std::cout << "value range: [" << min << "," << max << "]" << std::endl; - if (threshold == std::numeric_limits::max()) + if (threshold == std::numeric_limits::infinity()) ripser(std::move(dist), dim_max, threshold, modulus) .compute_barcodes(); else -- cgit v1.2.3 From b6b8b14a129107fc5bd70c999d551a6c9632881d Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sat, 26 Aug 2017 08:19:09 +0200 Subject: filter by persistence ratio --- ripser.cpp | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 74f9224..2f41b05 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -394,6 +394,7 @@ template class ripser { DistanceMatrix dist; index_t n, dim_max; value_t threshold; + float ratio; coefficient_t modulus; const binomial_coeff_table binomial_coeff; std::vector multiplicative_inverse; @@ -402,10 +403,10 @@ template class ripser { mutable std::vector::const_reverse_iterator> neighbor_end; public: - ripser(DistanceMatrix&& _dist, index_t _dim_max, value_t _threshold, coefficient_t _modulus) + 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), - modulus(_modulus), binomial_coeff(n, dim_max + 2), + 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 { @@ -781,7 +782,7 @@ public: found_persistence_pair: #ifdef PRINT_PERSISTENCE_PAIRS value_t death = get_diameter(pivot); - if (diameter != death) { + if ((float)death/diameter > ratio) { #ifdef INDICATE_PROGRESS std::cout << "\033[K"; #endif @@ -1003,6 +1004,9 @@ 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; @@ -1024,6 +1028,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") @@ -1080,11 +1089,11 @@ int main(int argc, char** argv) { << std::endl; if (threshold == std::numeric_limits::infinity()) - ripser(std::move(dist), dim_max, threshold, modulus) + ripser(std::move(dist), dim_max, threshold, ratio, modulus) .compute_barcodes(); else ripser(sparse_distance_matrix(std::move(dist), threshold), dim_max, - threshold, modulus) + threshold, ratio, modulus) .compute_barcodes(); } -- cgit v1.2.3 From 922acb40aec110347ac17da391ce5ae2b5da57db Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 24 Sep 2017 18:19:44 +0200 Subject: consolidated code for sparse and dense distance matrices --- ripser.cpp | 593 +++++++++++++++++++++++++++++-------------------------------- 1 file changed, 283 insertions(+), 310 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 2f41b05..1cb1d79 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -456,201 +456,12 @@ 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()); - } - - 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); - } - }; - - class simplex_sparse_coboundary_enumerator { - private: - const ripser& parent; - - index_t idx_below, idx_above, v, k, max_vertex_below; - const diameter_entry_t simplex; - const coefficient_t modulus; - const DistanceMatrix& dist; - const binomial_coeff_table& binomial_coeff; - - std::vector& vertices; - std::vector::const_reverse_iterator>& neighbor_it; - std::vector::const_reverse_iterator>& neighbor_end; - diameter_index_t x; - - public: - simplex_sparse_coboundary_enumerator(const diameter_entry_t _simplex, 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), - dist(parent.dist), binomial_coeff(parent.binomial_coeff), vertices(parent.vertices), - neighbor_it(parent.neighbor_it), neighbor_end(parent.neighbor_end) { - - neighbor_it.clear(); - neighbor_end.clear(); - vertices.clear(); - - parent.get_simplex_vertices(idx_below, _dim, parent.n, std::back_inserter(vertices)); - - for (auto v : vertices) { - neighbor_it.push_back(dist.neighbors[v].rbegin()); - neighbor_end.push_back(dist.neighbors[v].rend()); - } - } - - bool has_next(bool all_cofaces = true) { - for (auto &it0 = neighbor_it[0], &end0 = neighbor_end[0]; it0 != end0; ++it0) { - x = *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)) - if (++it == end) return false; - auto y = *it; - if (get_index(y) != get_index(x)) - goto continue_outer; - else - x = std::max(x, y); - } - return all_cofaces || - !(k > 0 && - parent.get_next_vertex(max_vertex_below, idx_below, k) > get_index(x)); - continue_outer:; - } - return false; - } - - 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)); - - 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); - } - }; - - 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(); - -#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 - } - - void assemble_sparse_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_sparse_coboundary_enumerator cofaces(simplex, 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 + class simplex_coboundary_enumerator; - std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), - greater_diameter_or_smaller_index()); -#ifdef INDICATE_PROGRESS - std::cout << "\033[K"; -#endif - } + void assemble_columns_to_reduce(std::vector& simplices, + std::vector& columns_to_reduce, + hash_map& pivot_column_index, + index_t dim); template void compute_pairs(std::vector& columns_to_reduce, @@ -829,11 +640,286 @@ public: #endif } - void compute_barcodes(); + std::vector get_edges(); + + void compute_barcodes() { + + std::vector columns_to_reduce; + + std::vector simplices; + + { + union_find dset(n); + std::vector& edges = simplices; + + 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) { + 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(simplices, columns_to_reduce, pivot_column_index, + dim + 1); + } + } + } }; -template <> void ripser::compute_barcodes(); -template <> void ripser::compute_barcodes(); + +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 <> +class ripser::simplex_coboundary_enumerator { +private: + const ripser& parent; + + index_t idx_below, idx_above, v, k, max_vertex_below; + 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; + diameter_index_t x; + +public: + simplex_coboundary_enumerator(const diameter_entry_t _simplex, 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), + dist(parent.dist), binomial_coeff(parent.binomial_coeff), vertices(parent.vertices), + neighbor_it(parent.neighbor_it), neighbor_end(parent.neighbor_end) { + + neighbor_it.clear(); + neighbor_end.clear(); + vertices.clear(); + + parent.get_simplex_vertices(idx_below, _dim, parent.n, std::back_inserter(vertices)); + + for (auto v : vertices) { + neighbor_it.push_back(dist.neighbors[v].rbegin()); + neighbor_end.push_back(dist.neighbors[v].rend()); + } + } + + bool has_next(bool all_cofaces = true) { + for (auto &it0 = neighbor_it[0], &end0 = neighbor_end[0]; it0 != end0; ++it0) { + x = *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)) + if (++it == end) return false; + auto y = *it; + if (get_index(y) != get_index(x)) + goto continue_outer; + else + x = std::max(x, y); + } + return all_cofaces || + !(k > 0 && + parent.get_next_vertex(max_vertex_below, idx_below, k) > get_index(x)); + continue_outer:; + } + return false; + } + + 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)); + + 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); + } +}; + + +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 <> std::vector ripser::get_edges() { + std::vector edges; + 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))); + } + 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(simplex, 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, @@ -1072,8 +1158,6 @@ 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()); - value_t min = std::numeric_limits::infinity(), max = -std::numeric_limits::infinity(); for (auto d: dist.distances) { @@ -1097,114 +1181,3 @@ int main(int argc, char** argv) { .compute_barcodes(); } - -template <> 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); - } - } -} - -template <> void ripser::compute_barcodes() { - - std::vector columns_to_reduce; - - std::vector simplices; - - { - union_find dset(n); - std::vector& edges = simplices; - 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))); - } - 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_sparse_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, - dim + 1); - } - } -} -- cgit v1.2.3 From 12bda394bdf6d2c610600b246b500908f88e4d10 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 24 Sep 2017 18:59:04 +0200 Subject: extracted dim 0 computation --- ripser.cpp | 80 +++++++++++++++++++++++++++++++------------------------------- 1 file changed, 40 insertions(+), 40 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 1cb1d79..688df6a 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -462,6 +462,42 @@ public: std::vector& columns_to_reduce, hash_map& pivot_column_index, index_t dim); + + void compute_dim_0_pairs(std::vector& edges, + std::vector& columns_to_reduce) { + 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) { + 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 + } template void compute_pairs(std::vector& columns_to_reduce, @@ -644,56 +680,20 @@ public: void compute_barcodes() { - std::vector columns_to_reduce; + std::vector simplices, columns_to_reduce; - std::vector simplices; - - { - union_find dset(n); - std::vector& edges = simplices; - - 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) { - 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 - } + 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); + dim); if (dim < dim_max) { assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, - dim + 1); + dim + 1); } } } -- cgit v1.2.3 From 740992e73b581c646f0f381789264f2504e76e8b Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 24 Sep 2017 19:16:08 +0200 Subject: use sparse distance matrix only when threshold is below maximum value --- ripser.cpp | 45 +++++++++++++++++++++++---------------------- 1 file changed, 23 insertions(+), 22 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 688df6a..5acc32e 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -1089,10 +1089,9 @@ int main(int argc, char** argv) { file_format format = DISTANCE_MATRIX; index_t dim_max = 1; - value_t threshold = std::numeric_limits::infinity(); - - float ratio = 1; + value_t threshold = std::numeric_limits::max(); + float ratio = 1; #ifdef USE_COEFFICIENTS coefficient_t modulus = 2; @@ -1156,28 +1155,30 @@ 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()); - } + 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; } - - std::cout << "value range: [" << min << "," << max << "]" - << std::endl; - if (threshold == std::numeric_limits::infinity()) - ripser(std::move(dist), dim_max, threshold, ratio, modulus) - .compute_barcodes(); - else + std::cout << "value range: [" << min << "," << max_finite << "]" << std::endl; + + if (threshold >= max) { + std::cout << "distance matrix with " << dist.size() << " points" << std::endl; + ripser(std::move(dist), dim_max, threshold, ratio, + modulus) + .compute_barcodes(); + } else { + std::cout << "sparse distance matrix with " << dist.size() << " points and " << num_edges + << " entries" << std::endl; + ripser(sparse_distance_matrix(std::move(dist), threshold), dim_max, threshold, ratio, modulus) - .compute_barcodes(); - + .compute_barcodes(); + } } -- cgit v1.2.3 From a962f06a93df50941c813d187730b849b8fd84d6 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 24 Sep 2017 19:16:35 +0200 Subject: pretty-print --- ripser.cpp | 190 +++++++++++++++++++++++++++++-------------------------------- 1 file changed, 91 insertions(+), 99 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 5acc32e..70a092a 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -146,7 +146,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; } @@ -249,10 +250,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()); @@ -403,7 +402,8 @@ template class ripser { mutable std::vector::const_reverse_iterator> neighbor_end; public: - ripser(DistanceMatrix&& _dist, index_t _dim_max, value_t _threshold, float _ratio, coefficient_t _modulus) + 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), ratio(_ratio), modulus(_modulus), binomial_coeff(n, dim_max + 2), @@ -459,29 +459,28 @@ 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); - + std::vector& columns_to_reduce, + hash_map& pivot_column_index, index_t dim); + void compute_dim_0_pairs(std::vector& edges, - std::vector& columns_to_reduce) { + std::vector& columns_to_reduce) { union_find dset(n); - + edges = get_edges(); - + std::sort(edges.rbegin(), edges.rend(), - greater_diameter_or_smaller_index()); - + 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) @@ -492,7 +491,7 @@ public: 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; @@ -629,7 +628,7 @@ public: found_persistence_pair: #ifdef PRINT_PERSISTENCE_PAIRS value_t death = get_diameter(pivot); - if ((float)death/diameter > ratio) { + if ((float)death / diameter > ratio) { #ifdef INDICATE_PROGRESS std::cout << "\033[K"; #endif @@ -677,31 +676,29 @@ public: } 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); - + dim); + if (dim < dim_max) { assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, - dim + 1); + dim + 1); } } } }; - -template <> -class ripser::simplex_coboundary_enumerator { +template <> class ripser::simplex_coboundary_enumerator { private: index_t idx_below, idx_above, v, k; std::vector vertices; @@ -709,16 +706,16 @@ private: 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) { + 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); @@ -729,53 +726,52 @@ 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; + (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; return diameter_entry_t(coface_diameter, coface_index, coface_coefficient); } }; -template <> -class ripser::simplex_coboundary_enumerator { +template <> class ripser::simplex_coboundary_enumerator { private: const ripser& parent; - + index_t idx_below, idx_above, v, k, max_vertex_below; 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; diameter_index_t x; - + public: simplex_coboundary_enumerator(const diameter_entry_t _simplex, 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), - dist(parent.dist), binomial_coeff(parent.binomial_coeff), vertices(parent.vertices), - neighbor_it(parent.neighbor_it), neighbor_end(parent.neighbor_end) { - + 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) { + neighbor_it.clear(); neighbor_end.clear(); vertices.clear(); - + parent.get_simplex_vertices(idx_below, _dim, parent.n, std::back_inserter(vertices)); - + for (auto v : vertices) { neighbor_it.push_back(dist.neighbors[v].rbegin()); neighbor_end.push_back(dist.neighbors[v].rend()); } } - + bool has_next(bool all_cofaces = true) { for (auto &it0 = neighbor_it[0], &end0 = neighbor_end[0]; it0 != end0; ++it0) { x = *it0; @@ -790,34 +786,33 @@ public: x = std::max(x, y); } return all_cofaces || - !(k > 0 && - parent.get_next_vertex(max_vertex_below, idx_below, k) > get_index(x)); + !(k > 0 && + parent.get_next_vertex(max_vertex_below, idx_below, k) > get_index(x)); continue_outer:; } return false; } - + 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)); - + coefficient_t coface_coefficient = - (k & 1 ? -1 + modulus : 1) * get_coefficient(simplex) % modulus; - + (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); + idx_above + binomial_coeff(get_index(x), k + 1) + idx_below, + coface_coefficient); } }; - template <> std::vector ripser::get_edges() { std::vector edges; for (index_t index = binomial_coeff(n, 2); index-- > 0;) { @@ -837,85 +832,82 @@ 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) { +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"; + << "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)); + 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"; + << "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"; + << "sorting " << num_simplices << " columns" << std::flush << "\r"; #endif - + std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), - greater_diameter_or_smaller_index()); + 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) { - +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"; + << "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(simplex, 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))); + 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"; + << "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()); + greater_diameter_or_smaller_index()); #ifdef INDICATE_PROGRESS std::cout << "\033[K"; #endif -- cgit v1.2.3 From 23a262f1b59eacd7a3ef01589101e64290cd7128 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 15 Feb 2018 11:00:13 +0100 Subject: threshold at enclosing radius --- ripser.cpp | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 70a092a..7858c55 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -785,9 +785,8 @@ 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)); + return all_cofaces || !(k > 0 && parent.get_next_vertex(max_vertex_below, idx_below, + k) > get_index(x)); continue_outer:; } return false; @@ -1151,6 +1150,15 @@ int main(int argc, char** argv) { max = -std::numeric_limits::infinity(), max_finite = max; int num_edges = 0; + 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); + } + + if (threshold == std::numeric_limits::max()) threshold = enclosing_radius; + for (auto d : dist.distances) { min = std::min(min, d); max = std::max(max, d); -- cgit v1.2.3