diff options
| -rw-r--r-- | ripser.cpp | 224 |
1 files changed, 151 insertions, 73 deletions
@@ -180,6 +180,7 @@ class diameter_index_t : public std::pair<value_t, index_t> { public: diameter_index_t() : std::pair<value_t, index_t>() {} diameter_index_t(std::pair<value_t, index_t> p) : std::pair<value_t, index_t>(p) {} + diameter_index_t(index_t i) : std::pair<value_t, index_t>(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<value_t, entry_t> { public: diameter_entry_t(std::pair<value_t, entry_t> p) : std::pair<value_t, entry_t>(p) {} - diameter_entry_t(entry_t e) : std::pair<value_t, entry_t>(0, e) {} - diameter_entry_t() : diameter_entry_t(0) {} diameter_entry_t(value_t _diameter, index_t _index, coefficient_t _coefficient) : std::pair<value_t, entry_t>(_diameter, make_entry(_index, _coefficient)) {} diameter_entry_t(diameter_index_t _diameter_index, coefficient_t _coefficient) @@ -211,42 +210,6 @@ template <typename Entry> struct greater_diameter_or_smaller_index { } }; -template <typename DistanceMatrix> class rips_filtration_comparator { -public: - const DistanceMatrix& dist; - const index_t dim; - -private: - mutable std::vector<index_t> 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_index_t(diameter(a), a), - diameter_index_t(diameter(b), b)); - } - - template <typename Entry> bool operator()(const Entry& a, const Entry& b) const { - return operator()(get_index(a), get_index(b)); - } -}; - template <class DistanceMatrix> class simplex_coboundary_enumerator { private: index_t idx_below, idx_above, v, k; @@ -317,6 +280,96 @@ public: size_t size() const { return rows.size(); } }; +class sparse_distance_matrix { +public: + std::vector<std::vector<diameter_index_t>> neighbors; + + template <typename DistanceMatrix> + 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(std::make_pair(mat(i, j), j)); + } + + size_t size() const { return neighbors.size(); } +}; + +template <class T> struct second_argument_greater { + bool operator()(const T& lhs, const T& rhs) const { return lhs.second > rhs.second; } +}; + +template <class T> struct second_argument_equal_to { + bool operator()(const T& lhs, const T& rhs) const { return lhs.second == rhs.second; } +}; + +template <> class simplex_coboundary_enumerator<const sparse_distance_matrix&> { +private: + 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<index_t> vertices; + + const diameter_entry_t simplex; + + std::vector<std::vector<diameter_index_t>::const_reverse_iterator> ii, ee; + + diameter_index_t x; + + const coefficient_t modulus; + +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) + : 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) { + 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; + 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<diameter_index_t>()); + } + return all_cofaces || + !(k > 0 && get_next_vertex(max_vertex_below, idx_below, k, binomial_coeff) > get_index(x)); + continue_outer:; + } + return false; + } + + 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); + 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 <> void compressed_distance_matrix<LOWER_TRIANGULAR>::init_rows() { value_t* pointer = &distances[0]; for (index_t i = 1; i < size(); ++i) { @@ -485,35 +538,41 @@ template <typename Heap> void push_entry(Heap& column, index_t i, coefficient_t column.push(std::make_pair(diameter, e)); } -template <typename Comparator> -void assemble_columns_to_reduce(std::vector<diameter_index_t>& columns_to_reduce, - hash_map<index_t, index_t>& 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(); +void assemble_columns_to_reduce(std::vector<diameter_index_t>& simplices, + std::vector<diameter_index_t>& columns_to_reduce, + hash_map<index_t, index_t>& pivot_column_index, + const sparse_distance_matrix& sparse_dist, 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" - << "assembling " << num_simplices << " columns" << std::flush << "\r"; + << "assembling 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 + columns_to_reduce.clear(); + + std::vector<diameter_index_t> next_simplices; + + for (diameter_index_t simplex : simplices) { + simplex_coboundary_enumerator<const sparse_distance_matrix&> 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 " << 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(), @@ -523,11 +582,10 @@ void assemble_columns_to_reduce(std::vector<diameter_index_t>& columns_to_reduce #endif } -template <typename DistanceMatrix, typename ComparatorCofaces, typename Comparator> +template <typename DistanceMatrix> void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce, hash_map<index_t, index_t>& pivot_column_index, index_t dim, index_t n, value_t threshold, coefficient_t modulus, const std::vector<coefficient_t>& multiplicative_inverse, const DistanceMatrix& dist, - const ComparatorCofaces& comp, const Comparator& comp_prev, const binomial_coeff_table& binomial_coeff) { #ifdef PRINT_PERSISTENCE_PAIRS @@ -909,8 +967,14 @@ int main(int argc, char** argv) { exit(-1); } + std::chrono::time_point<std::chrono::system_clock> 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::milliseconds>(std::chrono::system_clock::now() - start).count() / 1000. << " s.\n"; + index_t n = dist.size(); std::cout << "distance matrix with " << n << " points" << std::endl; @@ -918,6 +982,13 @@ 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::milliseconds>(std::chrono::system_clock::now() - start).count() / 1000. << " s.\n"; + + dim_max = std::min(dim_max, n - 2); binomial_coeff_table binomial_coeff(n, dim_max + 2); @@ -925,15 +996,20 @@ int main(int argc, char** argv) { std::vector<diameter_index_t> columns_to_reduce; - { + std::vector<diameter_index_t> simplices, &edges = simplices; + + + 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<diameter_index_t>()); + + { union_find dset(n); - std::vector<diameter_index_t> edges; - rips_filtration_comparator<decltype(dist)> 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<diameter_index_t>()); #ifdef PRINT_PERSISTENCE_PAIRS std::cout << "persistence intervals in dim 0:" << std::endl; @@ -962,17 +1038,19 @@ int main(int argc, char** argv) { } for (index_t dim = 1; dim <= dim_max; ++dim) { - rips_filtration_comparator<decltype(dist)> comp(dist, dim + 1, binomial_coeff); - rips_filtration_comparator<decltype(dist)> comp_prev(dist, dim, binomial_coeff); - hash_map<index_t, index_t> 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, - comp, comp_prev, binomial_coeff); + compute_pairs(columns_to_reduce, pivot_column_index, dim, n, threshold, modulus, multiplicative_inverse, + sparse_dist, 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, dim, n, + threshold, modulus, binomial_coeff); } } + + std::cout << "Computing Rips persistence in " << std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now() - start).count() / 1000. << " s.\n"; + + } |