diff options
Diffstat (limited to 'ripser.cpp')
| -rw-r--r-- | ripser.cpp | 161 |
1 files changed, 141 insertions, 20 deletions
@@ -94,30 +94,38 @@ 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 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; @@ -295,6 +303,111 @@ 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(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 <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 InputIteratorCollection, class Compare, class Equal> +class set_intersection_enumerator { +public: + InputIteratorCollection &first, &last; + 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; + 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]++; + } +}; + +class simplex_coboundary_enumerator_sparse { +private: + index_t idx_below, idx_above, v, k, w; + const binomial_coeff_table& binomial_coeff; + const sparse_distance_matrix& sparse_dist; + std::vector<index_t> vertices; + + std::vector<std::vector<diameter_index_t>::const_reverse_iterator> ii, ee; + set_intersection_enumerator<decltype(ii), second_argument_greater<diameter_index_t>, second_argument_equal_to<diameter_index_t>> 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_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()); + } + } + + bool has_next() { + return v_intersection.has_next(); + } + + std::pair<entry_t, index_t> 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); + --k; + } + + 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<LOWER_TRIANGULAR>::init_rows() { value_t* pointer = &distances[0]; for (index_t i = 1; i < size(); ++i) { @@ -469,6 +582,11 @@ void assemble_columns_to_reduce(std::vector<diameter_index_t>& 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 @@ -585,6 +703,9 @@ void compute_pairs(std::vector<diameter_index_t>& columns_to_reduce, hash_map<in coface_entries.clear(); simplex_coboundary_enumerator cofaces(get_index(simplex), dim, n, binomial_coeff); + + // iterate over cofaces using set intersection of neighbors + while (cofaces.has_next()) { auto coface_descriptor = cofaces.next(); entry_t coface = coface_descriptor.first; @@ -950,4 +1071,4 @@ int main(int argc, char** argv) { assemble_columns_to_reduce(columns_to_reduce, pivot_column_index, comp, dim, n, threshold, binomial_coeff); } } -}
\ No newline at end of file +} |