diff options
Diffstat (limited to 'ripser.cpp')
-rw-r--r-- | ripser.cpp | 629 |
1 files changed, 286 insertions, 343 deletions
@@ -47,35 +47,25 @@ template <class Key, class T> class hash_map : public std::unordered_map<Key, T> #endif typedef float value_t; -// typedef uint16_t value_t; - typedef int64_t index_t; typedef int16_t coefficient_t; class binomial_coeff_table { std::vector<std::vector<index_t>> B; - index_t n_max, k_max; - public: - binomial_coeff_table(index_t n, index_t k) { - n_max = n; - k_max = k; - - B.resize(n + 1); + binomial_coeff_table(index_t n, index_t k) : B(n + 1) { for (index_t i = 0; i <= n; i++) { B[i].resize(k + 1); - for (index_t j = 0; j <= std::min(i, k); j++) { + for (index_t j = 0; j <= std::min(i, k); j++) if (j == 0 || j == i) B[i][j] = 1; else B[i][j] = B[i - 1][j - 1] + B[i - 1][j]; - } } } index_t operator()(index_t n, index_t k) const { - assert(n <= n_max); - assert(k <= k_max); + assert(n < B.size() && k < B[n].size()); return B[n][k]; } }; @@ -97,43 +87,6 @@ std::vector<coefficient_t> multiplicative_inverse_vector(const coefficient_t m) return inverse; } -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; - 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); - return v; -} - -template <typename OutputIterator> -OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t v, - const binomial_coeff_table& binomial_coeff, OutputIterator out) { - --v; - for (index_t k = dim + 1; k > 0; --k) { - get_next_vertex(v, idx, k, binomial_coeff); - *out++ = v; - idx -= binomial_coeff(v, k); - } - return out; -} - -std::vector<index_t> vertices_of_simplex(const index_t simplex_index, const index_t dim, const index_t n, - const binomial_coeff_table& binomial_coeff) { - std::vector<index_t> vertices; - get_simplex_vertices(simplex_index, dim, n, binomial_coeff, std::back_inserter(vertices)); - return vertices; -} #ifdef USE_COEFFICIENTS struct __attribute__((packed)) entry_t { @@ -172,10 +125,6 @@ void set_coefficient(index_t& e, const coefficient_t c) {} const entry_t& get_entry(const entry_t& e) { return e; } -template <typename Entry> 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<value_t, index_t> { public: diameter_index_t() : std::pair<value_t, index_t>() {} @@ -210,47 +159,6 @@ template <typename Entry> struct greater_diameter_or_smaller_index { } }; -template <class DistanceMatrix> class simplex_coboundary_enumerator { -private: - index_t idx_below, idx_above, v, k; - std::vector<index_t> vertices; - const diameter_entry_t simplex; - const coefficient_t modulus; - const DistanceMatrix& dist; - const binomial_coeff_table& binomial_coeff; - -public: - 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)) { - idx_below -= binomial_coeff(v, k); - idx_above += binomial_coeff(v, k + 1); - - --v; - --k; - assert(k != -1); - } - return v != -1; - } - - 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; - return diameter_entry_t(coface_diameter, idx_above + binomial_coeff(v--, k + 1) + idx_below, - coface_coefficient); - } -}; - enum compressed_matrix_layout { LOWER_TRIANGULAR, UPPER_TRIANGULAR }; template <compressed_matrix_layout Layout> class compressed_distance_matrix { @@ -303,73 +211,6 @@ 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) { @@ -538,12 +379,135 @@ template <typename Heap> void push_entry(Heap& column, index_t i, coefficient_t column.push(std::make_pair(diameter, e)); } +class ripser { + index_t dim_max, n; + value_t threshold; + const binomial_coeff_table binomial_coeff; + std::vector<coefficient_t> multiplicative_inverse; + coefficient_t modulus; + compressed_lower_distance_matrix dist; + sparse_distance_matrix sparse_dist; + mutable std::vector<index_t> vertices; + +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)) {} + + 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<index_t> vertices; + + std::vector<std::vector<diameter_index_t>::const_reverse_iterator> ii, 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) { + 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; + 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 && 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); + } + }; + + + 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; + } + + template <typename OutputIterator> + 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; + } + +// value_t compute_diameter(const index_t index, index_t dim) const { +// value_t diam = -std::numeric_limits<value_t>::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; +// } + 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) { + hash_map<index_t, index_t>& pivot_column_index, index_t dim) { #ifdef INDICATE_PROGRESS std::cout << "\033[K" @@ -555,8 +519,7 @@ void assemble_columns_to_reduce(std::vector<diameter_index_t>& simplices, 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); + simplex_coboundary_enumerator cofaces(simplex, dim, *this); while (cofaces.has_next(false)) { auto coface = cofaces.next(); @@ -582,177 +545,232 @@ void assemble_columns_to_reduce(std::vector<diameter_index_t>& simplices, #endif } -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 binomial_coeff_table& binomial_coeff) { + void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce, hash_map<index_t, index_t>& pivot_column_index, + index_t dim) { #ifdef PRINT_PERSISTENCE_PAIRS - std::cout << "persistence intervals in dim " << dim << ":" << std::endl; + std::cout << "persistence intervals in dim " << dim << ":" << std::endl; #endif #ifdef ASSEMBLE_REDUCTION_MATRIX - compressed_sparse_matrix<diameter_entry_t> reduction_coefficients; + compressed_sparse_matrix<diameter_entry_t> reduction_coefficients; #else #ifdef USE_COEFFICIENTS - std::vector<diameter_entry_t> reduction_coefficients; + std::vector<diameter_entry_t> reduction_coefficients; #endif #endif - std::vector<diameter_entry_t> coface_entries; + std::vector<diameter_entry_t> coface_entries; - for (index_t i = 0; i < columns_to_reduce.size(); ++i) { - auto column_to_reduce = columns_to_reduce[i]; + for (index_t i = 0; i < columns_to_reduce.size(); ++i) { + auto column_to_reduce = columns_to_reduce[i]; #ifdef ASSEMBLE_REDUCTION_MATRIX - std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, greater_diameter_or_smaller_index<diameter_entry_t>> - reduction_column; + std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, + greater_diameter_or_smaller_index<diameter_entry_t>> + reduction_column; #endif - std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, - greater_diameter_or_smaller_index<diameter_entry_t>> - working_coboundary; + std::priority_queue<diameter_entry_t, std::vector<diameter_entry_t>, + greater_diameter_or_smaller_index<diameter_entry_t>> + working_coboundary; - value_t diameter = get_diameter(column_to_reduce); + value_t diameter = get_diameter(column_to_reduce); #ifdef INDICATE_PROGRESS - if ((i + 1) % 1000 == 0) - std::cout << "\033[K" - << "reducing column " << i + 1 << "/" << columns_to_reduce.size() << " (diameter " << diameter - << ")" << std::flush << "\r"; + if ((i + 1) % 1000 == 0) + std::cout << "\033[K" + << "reducing column " << i + 1 << "/" << columns_to_reduce.size() << " (diameter " << diameter + << ")" << std::flush << "\r"; #endif - index_t j = i; + index_t j = i; - // start with a dummy pivot entry with coefficient -1 in order to initialize - // working_coboundary with the coboundary of the simplex with index column_to_reduce - diameter_entry_t pivot(0, -1, -1 + modulus); + // start with a dummy pivot entry with coefficient -1 in order to initialize + // working_coboundary with the coboundary of the simplex with index column_to_reduce + diameter_entry_t pivot(0, -1, -1 + modulus); #ifdef ASSEMBLE_REDUCTION_MATRIX - // initialize reduction_coefficients as identity matrix - reduction_coefficients.append_column(); + // initialize reduction_coefficients as identity matrix + reduction_coefficients.append_column(); #endif #ifdef USE_COEFFICIENTS - reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, 1)); + reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, 1)); #endif - - bool might_be_apparent_pair = (i == j); - do { - const coefficient_t factor = modulus - get_coefficient(pivot); + bool might_be_apparent_pair = (i == j); + + do { + const coefficient_t factor = modulus - get_coefficient(pivot); #ifdef ASSEMBLE_REDUCTION_MATRIX #ifdef USE_COEFFICIENTS - auto coeffs_begin = reduction_coefficients.cbegin(j), coeffs_end = reduction_coefficients.cend(j); + auto coeffs_begin = reduction_coefficients.cbegin(j), coeffs_end = reduction_coefficients.cend(j); #else - std::vector<diameter_entry_t> coeffs(0); - coeffs.push_back(columns_to_reduce[j]); - for (auto it = reduction_coefficients.cbegin(j); it != reduction_coefficients.cend(j); ++it) coeffs.push_back(*it); - auto coeffs_begin = coeffs.begin(), coeffs_end = coeffs.end(); + std::vector<diameter_entry_t> coeffs; + coeffs.push_back(columns_to_reduce[j]); + for (auto it = reduction_coefficients.cbegin(j); it != reduction_coefficients.cend(j); ++it) + coeffs.push_back(*it); + auto coeffs_begin = coeffs.begin(), coeffs_end = coeffs.end(); #endif #else #ifdef USE_COEFFICIENTS - auto coeffs_begin = &reduction_coefficients[j], coeffs_end = &reduction_coefficients[j] + 1; + auto coeffs_begin = &reduction_coefficients[j], coeffs_end = &reduction_coefficients[j] + 1; #else - auto coeffs_begin = &columns_to_reduce[j], coeffs_end = &columns_to_reduce[j] + 1; + auto coeffs_begin = &columns_to_reduce[j], coeffs_end = &columns_to_reduce[j] + 1; #endif #endif - for (auto it = coeffs_begin; it != coeffs_end; ++it) { - diameter_entry_t simplex = *it; - set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); + for (auto it = coeffs_begin; it != coeffs_end; ++it) { + diameter_entry_t simplex = *it; + set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); #ifdef ASSEMBLE_REDUCTION_MATRIX - reduction_column.push(simplex); -#endif - - coface_entries.clear(); - simplex_coboundary_enumerator<decltype(dist)> cofaces(simplex, dim, n, modulus, dist, binomial_coeff); - while (cofaces.has_next()) { - diameter_entry_t coface = cofaces.next(); - if (get_diameter(coface) <= threshold) { - coface_entries.push_back(coface); - if (might_be_apparent_pair && (get_diameter(simplex) == get_diameter(coface))) { - if (pivot_column_index.find(get_index(coface)) == pivot_column_index.end()) { - pivot = coface; - goto found_persistence_pair; + reduction_column.push(simplex); +#endif + + coface_entries.clear(); + simplex_coboundary_enumerator cofaces(simplex, dim, *this); + while (cofaces.has_next()) { + diameter_entry_t coface = cofaces.next(); + if (get_diameter(coface) <= threshold) { + coface_entries.push_back(coface); + if (might_be_apparent_pair && (get_diameter(simplex) == get_diameter(coface))) { + if (pivot_column_index.find(get_index(coface)) == pivot_column_index.end()) { + pivot = coface; + goto found_persistence_pair; + } + might_be_apparent_pair = false; } - might_be_apparent_pair = false; } } + for (auto coface : coface_entries) working_coboundary.push(coface); } - for (auto coface : coface_entries) working_coboundary.push(coface); - } - pivot = get_pivot(working_coboundary, modulus); + pivot = get_pivot(working_coboundary, modulus); - if (get_index(pivot) != -1) { - auto pair = pivot_column_index.find(get_index(pivot)); + if (get_index(pivot) != -1) { + auto pair = pivot_column_index.find(get_index(pivot)); - if (pair != pivot_column_index.end()) { - j = pair->second; - continue; - } - } else { + if (pair != pivot_column_index.end()) { + j = pair->second; + continue; + } + } else { #ifdef PRINT_PERSISTENCE_PAIRS #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif - std::cout << " [" << diameter << ", )" << std::endl << std::flush; + std::cout << " [" << diameter << ", )" << std::endl << std::flush; #endif - break; - } + break; + } - found_persistence_pair: + found_persistence_pair: #ifdef PRINT_PERSISTENCE_PAIRS - value_t death = get_diameter(pivot); - if (diameter != death) { + value_t death = get_diameter(pivot); + if (diameter != death) { #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif - std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush; - } + std::cout << " [" << diameter << "," << death << ")" << std::endl << std::flush; + } #endif - pivot_column_index.insert(std::make_pair(get_index(pivot), i)); + pivot_column_index.insert(std::make_pair(get_index(pivot), i)); #ifdef USE_COEFFICIENTS - const coefficient_t inverse = multiplicative_inverse[get_coefficient(pivot)]; + const coefficient_t inverse = multiplicative_inverse[get_coefficient(pivot)]; #endif #ifdef ASSEMBLE_REDUCTION_MATRIX - // replace current column of reduction_coefficients (with a single diagonal 1 entry) - // by reduction_column (possibly with a different entry on the diagonal) +// replace current column of reduction_coefficients (with a single diagonal 1 entry) +// by reduction_column (possibly with a different entry on the diagonal) #ifdef USE_COEFFICIENTS - reduction_coefficients.pop_back(); + reduction_coefficients.pop_back(); #else - pop_pivot(reduction_column, modulus); + pop_pivot(reduction_column, modulus); #endif - - while (true) { - diameter_entry_t e = pop_pivot(reduction_column, modulus); - if (get_index(e) == -1) break; + + while (true) { + diameter_entry_t e = pop_pivot(reduction_column, modulus); + if (get_index(e) == -1) break; #ifdef USE_COEFFICIENTS - set_coefficient(e, inverse * get_coefficient(e) % modulus); - assert(get_coefficient(e) > 0); + set_coefficient(e, inverse * get_coefficient(e) % modulus); + assert(get_coefficient(e) > 0); #endif - reduction_coefficients.push_back(e); - } + reduction_coefficients.push_back(e); + } #else #ifdef USE_COEFFICIENTS - reduction_coefficients.pop_back(); - reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, inverse)); + reduction_coefficients.pop_back(); + reduction_coefficients.push_back(diameter_entry_t(column_to_reduce, inverse)); #endif #endif - break; - } while (true); - } + break; + } while (true); + } #ifdef INDICATE_PROGRESS - std::cout << "\033[K"; + std::cout << "\033[K"; #endif -} + } + + void compute_barcodes() { + + std::vector<diameter_index_t> columns_to_reduce; + + std::vector<diameter_index_t> 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)); + } + + std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index<diameter_index_t>()); + + { + union_find dset(n); + +#ifdef PRINT_PERSISTENCE_PAIRS + std::cout << "persistence intervals in dim 0:" << std::endl; +#endif + + std::vector<index_t> 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<index_t, index_t> 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); + } + } + + } +}; enum file_format { LOWER_DISTANCE_MATRIX, UPPER_DISTANCE_MATRIX, DISTANCE_MATRIX, POINT_CLOUD, DIPHA, RIPSER }; @@ -967,90 +985,15 @@ 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; + std::cout << "distance matrix with " << dist.size() << " points" << std::endl; auto value_range = std::minmax_element(dist.distances.begin(), dist.distances.end()); std::cout << "value range: [" << *value_range.first << "," << *value_range.second << "]" << std::endl; - - 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); - std::vector<coefficient_t> multiplicative_inverse(multiplicative_inverse_vector(modulus)); - - 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); - -#ifdef PRINT_PERSISTENCE_PAIRS - std::cout << "persistence intervals in dim 0:" << std::endl; -#endif - - std::vector<index_t> vertices_of_edge(2); - for (auto e : edges) { - vertices_of_edge.clear(); - get_simplex_vertices(get_index(e), 1, n, binomial_coeff, 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<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, - sparse_dist, binomial_coeff); - - if (dim < dim_max) { - assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index, sparse_dist, dim, n, - threshold, modulus, binomial_coeff); - } - } + sparse_distance_matrix sparse_dist(dist, threshold); - std::cout << "Computing Rips persistence in " << std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now() - start).count() / 1000. << " s.\n"; - + ripser(std::move(dist), std::move(sparse_dist), dim_max, threshold, modulus).compute_barcodes(); } |