72#ifdef GUDHI_INDICATE_PROGRESS
76#include <boost/range/iterator_range_core.hpp>
77#include <boost/version.hpp>
78#if BOOST_VERSION >= 108100
79#include <boost/unordered/unordered_flat_map.hpp>
81#include <boost/unordered_map.hpp>
84#ifdef GUDHI_RIPSER_USE_BOOST_HEAP
85#include <boost/heap/d_ary_heap.hpp>
88#include <gudhi/uint128.h>
89#include <gudhi/Debug_utils.h>
92namespace Gudhi::ripser {
94#define GUDHI_assert(X) GUDHI_CHECK(X, std::logic_error(""))
96#ifdef GUDHI_RIPSER_USE_BOOST_HEAP
97template <
class T,
class,
class C>
98using Heap = boost::heap::d_ary_heap<T, boost::heap::arity<8>, boost::heap::compare<C>>;
100template <
class T,
class V,
class C>
101struct Heap : std::priority_queue<T, V, C> {
102 typedef std::priority_queue<T, V, C> Base;
104 void clear() { this->c.clear(); }
108template<
class vertex_t_>
111 typedef vertex_t_ vertex_t;
113 std::vector<vertex_t> parent;
114 std::vector<uint8_t> rank;
116 Union_find(
const vertex_t n) : parent(n), rank(n, 0) {
117 for (vertex_t i = 0; i < n; ++i) parent[i] = i;
121 vertex_t find(vertex_t x) {
124 while ((z = parent[y]) != y) y = z;
125 while ((z = parent[x]) != y) {
132 vertex_t y = parent[x];
133 vertex_t z = parent[y];
145 void link(vertex_t x, vertex_t y) {
147 if ((x = find(x)) == (y = find(y)))
return;
148 if (rank[x] > rank[y])
152 if (rank[x] == rank[y]) ++rank[y];
157template<
class coefficient_t>
158bool is_prime(
const coefficient_t n) {
159 if (!(n & 1) || n < 2)
return n == 2;
160 for (coefficient_t p = 3; p * p <= n; p += 2)
161 if (!(n % p))
return false;
166template<
class coefficient_t>
167coefficient_t normalize(
const coefficient_t n,
const coefficient_t modulus) {
168 return n > modulus / 2 ? n - modulus : n;
171template<
class coefficient_storage_t,
class coefficient_t>
172std::vector<coefficient_storage_t> multiplicative_inverse_vector(
const coefficient_t m) {
173 std::vector<coefficient_storage_t> inverse(m);
175 throw std::domain_error(
"Modulus must be a prime number");
176 if ((m - 1) != (coefficient_storage_t)(m - 1))
177 throw std::overflow_error(
"Modulus is too large");
182 for (coefficient_t a = 2; a < m; ++a) inverse[a] = m - (inverse[m % a] * (m / a)) % m;
186#ifdef GUDHI_INDICATE_PROGRESS
187constexpr std::chrono::milliseconds time_step(40);
188constexpr const char* clear_line=
"\r\033[K";
208template <
class Params>
209struct Full_distance_matrix {
211 typedef Tag_dense Category;
212 typedef typename Params::vertex_t vertex_t;
213 typedef typename Params::value_t value_t;
214 std::vector<value_t> distances;
221 template <
typename DistanceMatrix>
222 Full_distance_matrix(
const DistanceMatrix& mat)
223 : distances(static_cast<
std::size_t>(mat.size()) * mat.size()), n(mat.size()) {
224 for (vertex_t i = 0; i < size(); ++i)
225 for (vertex_t j = 0; j < size(); ++j)
226 distances[i * n + j] = mat(i, j);
230 value_t operator()(
const vertex_t j,
const vertex_t i)
const {
231 return distances[i * n + j];
233 vertex_t size()
const {
return n; }
236enum Compressed_matrix_layout { LOWER_TRIANGULAR, UPPER_TRIANGULAR };
238template <
class Params, Compressed_matrix_layout Layout>
239struct Compressed_distance_matrix {
241 typedef Tag_dense Category;
242 typedef typename Params::vertex_t vertex_t;
243 typedef typename Params::value_t value_t;
244 std::vector<value_t> distances;
246 std::vector<value_t*> rows;
249 Compressed_distance_matrix(std::vector<value_t>&& _distances)
250 : distances(
std::move(_distances)), rows((1 +
std::sqrt(1 + 8 * distances.size())) / 2) {
251 GUDHI_assert(distances.size() == (std::size_t)size() * (size() - 1) / 2);
255 template <
typename DistanceMatrix>
256 Compressed_distance_matrix(
const DistanceMatrix& mat)
257 : distances(static_cast<
std::size_t>(mat.size()) * (mat.size() - 1) / 2), rows(mat.size()) {
260 for (vertex_t i = 1; i < size(); ++i)
261 for (vertex_t j = 0; j < i; ++j) rows[i][j] = mat(i, j);
264 value_t operator()(
const vertex_t i,
const vertex_t j)
const {
265 if (i == j)
return 0;
266 if ((Layout == LOWER_TRIANGULAR) ? (i < j) : (i > j))
271 vertex_t size()
const {
return rows.size(); }
274 if constexpr (Layout == LOWER_TRIANGULAR) {
275 value_t* pointer = &distances[0];
276 for (vertex_t i = 1; i < size(); ++i) {
281 value_t* pointer = &distances[0] - 1;
282 for (vertex_t i = 0; i < size() - 1; ++i) {
284 pointer += size() - i - 2;
290template <
class Params>
291struct Sparse_distance_matrix {
293 typedef Tag_sparse Category;
294 static constexpr bool is_sparse =
true;
295 typedef typename Params::vertex_t vertex_t;
296 typedef typename Params::value_t value_t;
297 struct vertex_diameter_t {
298 vertex_diameter_t() =
default;
299 vertex_diameter_t(vertex_t i_, value_t d_) : i(i_), d(d_) {}
302 vertex_diameter_t(vertex_diameter_t
const&o)noexcept :i(o.i),d(o.d){}
303 vertex_diameter_t&operator=(vertex_diameter_t
const&o)
noexcept{i=o.i;d=o.d;
return*
this;}
305 vertex_t i; value_t d;
306 friend vertex_t get_vertex(
const vertex_diameter_t& i) {
return i.i; }
307 friend value_t get_diameter(
const vertex_diameter_t& i) {
return i.d; }
308 friend bool operator<(vertex_diameter_t
const& a, vertex_diameter_t
const& b) {
309 if (a.i < b.i)
return true;
310 if (a.i > b.i)
return false;
315 std::vector<std::vector<vertex_diameter_t>> neighbors;
316 std::size_t num_edges;
319#ifdef GUDHI_RIPSER_USE_HASHMAP_FOR_SPARSE_DIST_MAT
320#if BOOST_VERSION >= 108100
321 boost::unordered_flat_map<std::pair<vertex_t,vertex_t>,value_t> m;
323 boost::unordered_map<std::pair<vertex_t,vertex_t>,value_t> m;
329 Sparse_distance_matrix(std::vector<std::vector<vertex_diameter_t>>&& _neighbors,
330 std::size_t _num_edges = 0)
331 : neighbors(
std::move(_neighbors)), num_edges(_num_edges) {init();}
333 template <
typename DistanceMatrix>
334 Sparse_distance_matrix(
const DistanceMatrix& mat,
const value_t threshold)
335 : neighbors(mat.size()), num_edges(0) {
337 for (vertex_t i = 0; i < size(); ++i)
338 for (vertex_t j = 0; j < size(); ++j)
341 if (d <= threshold) {
343 neighbors[i].emplace_back(j, d);
349 value_t operator()(
const vertex_t i,
const vertex_t j)
const {
350#ifdef GUDHI_RIPSER_USE_HASHMAP_FOR_SPARSE_DIST_MAT
352 return m.at(std::minmax(i,j));
358 std::lower_bound(neighbors[i].begin(), neighbors[i].end(), vertex_diameter_t{j, 0});
359 return (neighbor != neighbors[i].end() && get_vertex(*neighbor) == j)
360 ? get_diameter(*neighbor)
361 :
std::numeric_limits<value_t>::infinity();
364 vertex_t size()
const {
return neighbors.size(); }
367#ifdef GUDHI_RIPSER_USE_HASHMAP_FOR_SPARSE_DIST_MAT
368 for(vertex_t i=0; i<size(); ++i){
369 for(
auto n:neighbors[i]){
370 m[std::minmax(i,get_vertex(n))]=get_diameter(n);
378template <
class Params>
379struct Euclidean_distance_matrix {
381 typedef Tag_other Category;
382 typedef typename Params::vertex_t vertex_t;
383 typedef typename Params::value_t value_t;
385 Euclidean_distance_matrix(std::vector<std::vector<value_t>>&& _points)
386 : points(
std::move(_points)) {
387 for (
auto p : points) { GUDHI_assert(p.size() == points.front().size()); }
390 value_t operator()(
const vertex_t i,
const vertex_t j)
const {
391 GUDHI_assert((std::size_t)i < points.size());
392 GUDHI_assert((std::size_t)j < points.size());
393 return std::sqrt(std::inner_product(
394 points[i].begin(), points[i].end(), points[j].begin(), value_t(), std::plus<value_t>(),
395 [](value_t u, value_t v) {
return (u - v) * (u - v); }));
398 vertex_t size()
const {
return points.size(); }
400 std::vector<std::vector<value_t>> points;
404template <
class DistanceMatrix,
class=std::
bool_constant<true>>
struct Is_sparse_impl : std::bool_constant<false> {};
405template <
class DistanceMatrix>
struct Is_sparse_impl<DistanceMatrix,
std::bool_constant<DistanceMatrix::is_sparse>> : std::bool_constant<true> {};
406template <
class DistanceMatrix>
constexpr bool is_sparse () {
return Is_sparse_impl<DistanceMatrix>::value; }
408template <
typename ValueType>
class Compressed_sparse_matrix_ {
409 std::vector<size_t> bounds;
410 std::vector<ValueType> entries;
412 typedef typename std::vector<ValueType>::const_iterator iterator;
413 typedef boost::iterator_range<iterator> iterator_pair;
416 iterator_pair subrange(
const size_t index)
const {
417 return {entries.begin() + (index == 0 ? 0 : bounds[index - 1]),
418 entries.begin() + bounds[index]};
422 void append_column() { bounds.push_back(entries.size()); }
424 void push_back(
const ValueType e) {
425 GUDHI_assert(0 < bounds.size());
426 entries.push_back(e);
444template <
class vertex_t>
445constexpr int log2up(vertex_t n) {
448 while(n>0) { n>>=1; ++k; }
452template <
class Params>
455 typedef typename Params::dimension_t dimension_t;
456 typedef typename Params::vertex_t vertex_t;
457 typedef typename Params::simplex_t simplex_t;
460 std::vector<std::vector<simplex_t>> B;
464 Cns_encoding(vertex_t n, dimension_t k) : B(k + 1,
std::vector<simplex_t>(n + 1, 0)) {
465 static_assert(std::numeric_limits<simplex_t>::radix == 2);
466 const int available_bits = std::numeric_limits<simplex_t>::digits;
467 simplex_t max_simplex_index = 0;
468 for (vertex_t i = 0; i <= n; ++i) {
470 for (dimension_t j = 1; (vertex_t)j < std::min<vertex_t>(i, k + 1); ++j)
471 B[j][i] = B[j - 1][i - 1] + B[j][i - 1];
472 if (i <= k) B[i][i] = 1;
473 vertex_t mi = std::min<vertex_t>(i >> 1, (vertex_t)k);
474 max_simplex_index = B[mi][i];
475 if (i > 1 && max_simplex_index < B[mi][i-1]) {
476 throw std::overflow_error(
"cannot encode all simplices of dimension " + std::to_string(k) +
" with " + std::to_string(n) +
" vertices using only " + std::to_string(available_bits) +
" bits");
479 extra_bits = available_bits - log2up(max_simplex_index + 1);
482 simplex_t operator()(vertex_t n, dimension_t k)
const {
483 GUDHI_assert(n >= k - 1);
488 vertex_t get_max_vertex(
const simplex_t idx,
const dimension_t k,
const vertex_t n)
const {
489 return get_max(n, k - 1, [&](vertex_t w) ->
bool {
return (*
this)(w, k) <= idx; });
492 int num_extra_bits()
const {
return extra_bits; }
495 template <
class Predicate>
496 static vertex_t get_max(vertex_t top,
const vertex_t bottom,
const Predicate pred) {
498 vertex_t count = top - bottom;
500 vertex_t step = count >> 1, mid = top - step;
512template <
class Params>
513class Bitfield_encoding {
515 typedef typename Params::dimension_t dimension_t;
516 typedef typename Params::vertex_t vertex_t;
517 typedef typename Params::simplex_t simplex_t;
524 Bitfield_encoding(vertex_t n, dimension_t k) : bits_per_vertex(log2up(n)) {
525 static_assert(std::numeric_limits<simplex_t>::radix == 2);
526 const int available_bits = std::numeric_limits<simplex_t>::digits;
527 extra_bits = available_bits - bits_per_vertex * k;
529 throw std::overflow_error(
"cannot encode all simplices of dimension " + std::to_string(k - 1) +
" with " + std::to_string(n) +
" vertices using only " + std::to_string(available_bits) +
" bits");
533 simplex_t operator()(vertex_t n, dimension_t k)
const {
538 return (simplex_t)(n - k) << (bits_per_vertex * k);
540 return (simplex_t)n << (bits_per_vertex * k);
544 vertex_t get_max_vertex(
const simplex_t idx, dimension_t k,
const vertex_t)
const {
548 return static_cast<vertex_t
>(idx >> (bits_per_vertex * k)) + k;
550 return static_cast<vertex_t
>(idx >> (bits_per_vertex * k));
554 int num_extra_bits()
const {
return extra_bits; }
557template <
typename DistanceMatrix,
typename SimplexEncoding,
typename Params>
struct Rips_filtration {
558 using size_t =
typename Params::size_t;
559 using vertex_t =
typename SimplexEncoding::vertex_t;
561 using simplex_t =
typename SimplexEncoding::simplex_t;
562 using dimension_t =
typename SimplexEncoding::dimension_t;
563 using value_t =
typename DistanceMatrix::value_t;
564 using coefficient_storage_t =
typename Params::coefficient_storage_t;
565 using coefficient_t =
typename Params::coefficient_t;
566 static constexpr bool use_coefficients = Params::use_coefficients;
569 struct entry_with_coeff_t {
571 entry_with_coeff_t(simplex_t _index, coefficient_t _coefficient,
int bits_for_coeff)
572 : content((_index << bits_for_coeff) | (_coefficient - 1)) { GUDHI_assert(_coefficient != 0); }
573 entry_with_coeff_t() {}
575 friend const entry_with_coeff_t& get_entry(
const entry_with_coeff_t& e) {
return e; }
577 simplex_t get_index(
const entry_with_coeff_t& e)
const {
return e.content >> num_bits_for_coeff(); }
578 coefficient_t get_coefficient(
const entry_with_coeff_t& e)
const {
return static_cast<coefficient_t
>(e.content & (((simplex_t)1 << num_bits_for_coeff()) - 1)) + 1; }
579 void set_coefficient(entry_with_coeff_t& e,
const coefficient_t c)
const { GUDHI_assert(c!=0); e.content = (e.content & ((simplex_t)(-1) << num_bits_for_coeff())) | (c - 1); }
582 struct entry_plain_t {
584 entry_plain_t(simplex_t _index, coefficient_t,
int) : index(_index) {}
586 friend const entry_plain_t& get_entry(
const entry_plain_t& e) {
return e; }
588 static simplex_t get_index(
const entry_plain_t& i) {
return i.index; }
589 static coefficient_t get_coefficient(
const entry_plain_t& i) {
return 1; }
590 static void set_coefficient(entry_plain_t& e,
const coefficient_t c) { GUDHI_assert(c==1); }
592 typedef std::conditional_t<use_coefficients, entry_with_coeff_t, entry_plain_t> entry_t;
593 entry_t make_entry(simplex_t i, coefficient_t c)
const {
return entry_t(i, c, num_bits_for_coeff()); }
595 static_assert(
sizeof(entry_t) ==
sizeof(simplex_t),
"size of entry_t is not the same as simplex_t");
599 Entry_hash(Rips_filtration
const& filt) : filtp(&filt) {}
600 Rips_filtration
const* filtp;
601 std::size_t operator()(
const entry_t& e)
const {
return boost::hash<simplex_t>()(filtp->get_index(e)); }
605 Equal_index(Rips_filtration
const& filt) : filtp(&filt) {}
606 Rips_filtration
const* filtp;
607 bool operator()(
const entry_t& e,
const entry_t& f)
const {
608 return filtp->get_index(e) == filtp->get_index(f);
612 struct diameter_simplex_t {
615 friend value_t get_diameter(
const diameter_simplex_t& i) {
return i.diameter; }
617 static simplex_t get_index(
const diameter_simplex_t& i) {
return i.index; }
619 struct diameter_entry_t : std::pair<value_t, entry_t> {
620 using std::pair<value_t, entry_t>::pair;
621 friend const entry_t& get_entry(
const diameter_entry_t& p) {
return p.second; }
622 friend entry_t& get_entry(diameter_entry_t& p) {
return p.second; }
623 friend value_t get_diameter(
const diameter_entry_t& p) {
return p.first; }
625 simplex_t get_index(
const diameter_entry_t& p)
const {
return get_index(get_entry(p)); }
626 coefficient_t get_coefficient(
const diameter_entry_t& p)
const {
627 return get_coefficient(get_entry(p));
629 void set_coefficient(diameter_entry_t& p,
const coefficient_t c)
const {
630 set_coefficient(get_entry(p), c);
632 diameter_entry_t make_diameter_entry(value_t diameter, simplex_t index, coefficient_t coefficient)
const {
633 return diameter_entry_t(diameter, make_entry(index, coefficient));
635 diameter_entry_t make_diameter_entry(
const diameter_simplex_t& diameter_index, coefficient_t coefficient)
const {
636 return diameter_entry_t(get_diameter(diameter_index), make_entry(get_index(diameter_index), coefficient));
639 template <
typename Entry>
struct Greater_diameter_or_smaller_index {
640 Greater_diameter_or_smaller_index(Rips_filtration
const& filt) : filtp(&filt) {}
641 Rips_filtration
const* filtp;
642 bool operator()(
const Entry& a,
const Entry& b)
const {
643 return (get_diameter(a) > get_diameter(b)) ||
644 ((get_diameter(a) == get_diameter(b)) && (filtp->get_index(a) < filtp->get_index(b)));
648 const DistanceMatrix dist;
650 const dimension_t dim_max;
651 const value_t threshold;
652 const coefficient_t modulus;
653 const SimplexEncoding simplex_encoding;
654 mutable std::vector<vertex_t> vertices;
657 Rips_filtration(DistanceMatrix&& _dist, dimension_t _dim_max, value_t _threshold, coefficient_t _modulus)
658 : dist(
std::move(_dist)), n(dist.size()),
659 dim_max(
std::min<vertex_t>(_dim_max, dist.size() - 2)), threshold(_threshold),
660 modulus(_modulus), simplex_encoding(n, dim_max + 2), bits_for_coeff(log2up(modulus-1)) {
662 if (use_coefficients && simplex_encoding.num_extra_bits() < num_bits_for_coeff())
664 throw std::overflow_error(
"Not enough spare bits in the simplex encoding to store a coefficient");
667 vertex_t num_vertices()
const {
return n; }
668 int num_bits_for_coeff()
const {
return bits_for_coeff; }
670 simplex_t get_edge_index(
const vertex_t i,
const vertex_t j)
const {
671 return simplex_encoding(i, 2) + j;
674 template <
typename OutputIterator>
675 OutputIterator get_simplex_vertices(simplex_t idx,
const dimension_t dim, vertex_t n,
676 OutputIterator out)
const {
678 for (dimension_t k = dim + 1; k > 1; --k) {
679 n = simplex_encoding.get_max_vertex(idx, k, n);
681 idx -= simplex_encoding(n, k);
683 *out =
static_cast<vertex_t
>(idx);
687 value_t compute_diameter(
const simplex_t index,
const dimension_t dim)
const {
688 value_t diam = -std::numeric_limits<value_t>::infinity();
690 vertices.resize(dim + 1);
691 get_simplex_vertices(index, dim, dist.size(), vertices.rbegin());
693 for (dimension_t i = 0; i <= dim; ++i)
694 for (dimension_t j = 0; j < i; ++j) {
695 diam = std::max(diam, dist(vertices[i], vertices[j]));
700 std::vector<diameter_simplex_t> get_edges() {
701 if constexpr (!std::is_same_v<typename DistanceMatrix::Category, Tag_sparse>) {
703 std::vector<diameter_simplex_t> edges;
704 for (vertex_t i = 0; i < n; ++i) {
705 for (vertex_t j = 0; j < i; ++j) {
706 value_t length = dist(i, j);
707 if (length <= threshold) edges.push_back({length, get_edge_index(i, j)});
712 std::vector<diameter_simplex_t> edges;
713 for (vertex_t i = 0; i < n; ++i)
714 for (
auto n : dist.neighbors[i]) {
715 vertex_t j = get_vertex(n);
716 if (i > j) edges.push_back({get_diameter(n), get_edge_index(i, j)});
723 template<
class DistanceMatrix2,
class=
typename DistanceMatrix2::Category>
class Simplex_coboundary_enumerator_ {
724 simplex_t idx_below, idx_above;
727 std::vector<vertex_t> vertices;
728 diameter_entry_t simplex;
729 const DistanceMatrix2& dist;
730 const SimplexEncoding& simplex_encoding;
731 const Rips_filtration& parent;
735 Simplex_coboundary_enumerator_(
const Rips_filtration& _parent) : dist(_parent.dist),
736 simplex_encoding(_parent.simplex_encoding), parent(_parent) {}
738 void set_simplex(
const diameter_entry_t _simplex,
const dimension_t _dim) {
739 idx_below = parent.get_index(_simplex);
744 vertices.resize(_dim + 1);
745 parent.get_simplex_vertices(parent.get_index(_simplex), _dim, dist.size(), vertices.rbegin());
748 bool has_next(
bool all_cofacets =
true) {
749 return (j >= k && (all_cofacets || simplex_encoding(j, k) > idx_below));
752 std::optional<diameter_entry_t> next_raw(
bool all_cofacets =
true) {
753 if (!has_next(all_cofacets))
return std::nullopt;
755 while (simplex_encoding(j, k) <= idx_below) {
756 idx_below -= simplex_encoding(j, k);
757 idx_above += simplex_encoding(j, k + 1);
760 GUDHI_assert(k != -1);
762 value_t cofacet_diameter = get_diameter(simplex);
764 for (vertex_t i : vertices) cofacet_diameter =
std::max(cofacet_diameter, dist(j, i));
765 simplex_t cofacet_index = idx_above + simplex_encoding(j--, k + 1) + idx_below;
766 coefficient_t cofacet_coefficient = parent.get_coefficient(simplex);
767 if (k & 1) cofacet_coefficient = parent.modulus - cofacet_coefficient;
768 return parent.make_diameter_entry(cofacet_diameter, cofacet_index, cofacet_coefficient);
770 std::optional<diameter_entry_t> next(
bool all_cofacets =
true) {
772 std::optional<diameter_entry_t> res = next_raw(all_cofacets);
773 if(!res || get_diameter(*res) <= parent.threshold)
return res;
778 template <
class DistanceMatrix2>
class Simplex_coboundary_enumerator_<DistanceMatrix2, Tag_sparse> {
779 typedef typename DistanceMatrix2::vertex_diameter_t vertex_diameter_t;
780 simplex_t idx_below, idx_above;
782 std::vector<vertex_t> vertices;
783 diameter_entry_t simplex;
784 const DistanceMatrix2& dist;
785 const SimplexEncoding& simplex_encoding;
786 std::vector<typename std::vector<vertex_diameter_t>::const_reverse_iterator> neighbor_it;
787 std::vector<typename std::vector<vertex_diameter_t>::const_reverse_iterator> neighbor_end;
788 vertex_diameter_t neighbor;
789 const Rips_filtration& parent;
792 Simplex_coboundary_enumerator_(
const Rips_filtration& _parent)
793 : dist(_parent.dist),
794 simplex_encoding(_parent.simplex_encoding), parent(_parent) {}
796 void set_simplex(
const diameter_entry_t _simplex,
const dimension_t _dim) {
797 idx_below = parent.get_index(_simplex);
801 vertices.resize(_dim + 1);
802 parent.get_simplex_vertices(idx_below, _dim, dist.size(), vertices.rbegin());
804 neighbor_it.resize(_dim + 1);
805 neighbor_end.resize(_dim + 1);
806 for (dimension_t i = 0; i <= _dim; ++i) {
807 auto v = vertices[i];
808 neighbor_it[i] = dist.neighbors[v].rbegin();
809 neighbor_end[i] = dist.neighbors[v].rend();
813 bool has_next(
bool all_cofacets =
true) {
814 for (
auto &it0 = neighbor_it[0], &end0 = neighbor_end[0]; it0 != end0; ++it0) {
816 for (
size_t idx = 1; idx < neighbor_it.size(); ++idx) {
817 auto &it = neighbor_it[idx], end = neighbor_end[idx];
818 while (get_vertex(*it) > get_vertex(neighbor))
819 if (++it == end)
return false;
820 if (get_vertex(*it) != get_vertex(neighbor))
823 neighbor = std::max(neighbor, *it);
825 while (k > 0 && vertices[k - 1] > get_vertex(neighbor)) {
826 if (!all_cofacets)
return false;
827 idx_below -= simplex_encoding(vertices[k - 1], k);
828 idx_above += simplex_encoding(vertices[k - 1], k + 1);
837 std::optional<diameter_entry_t> next_raw(
bool all_cofacets =
true) {
838 if (!has_next(all_cofacets))
return std::nullopt;
840 value_t cofacet_diameter = std::max(get_diameter(simplex), get_diameter(neighbor));
841 simplex_t cofacet_index = idx_above + simplex_encoding(get_vertex(neighbor), k + 1) + idx_below;
842 const coefficient_t modulus = parent.modulus;
843 coefficient_t cofacet_coefficient =
844 (k & 1 ? modulus - 1 : 1) * parent.get_coefficient(simplex) % modulus;
845 return parent.make_diameter_entry(cofacet_diameter, cofacet_index, cofacet_coefficient);
847 std::optional<diameter_entry_t> next(
bool all_cofacets =
true) {
848 return next_raw(all_cofacets);
852 typedef Simplex_coboundary_enumerator_<DistanceMatrix> Simplex_coboundary_enumerator;
854 class Simplex_boundary_enumerator {
856 simplex_t idx_below, idx_above;
859 diameter_entry_t simplex;
861 const SimplexEncoding& simplex_encoding;
862 const Rips_filtration& parent;
865 Simplex_boundary_enumerator(
const Rips_filtration& _parent)
866 : simplex_encoding(_parent.simplex_encoding), parent(_parent) {}
868 void set_simplex(
const diameter_entry_t _simplex,
const dimension_t _dim) {
869 idx_below = parent.get_index(_simplex);
877 bool has_next() {
return (k >= 0); }
879 std::optional<diameter_entry_t> next() {
880 if (!has_next())
return std::nullopt;
881 j = simplex_encoding.get_max_vertex(idx_below, k + 1, j);
883 simplex_t face_index = idx_above - simplex_encoding(j, k + 1) + idx_below;
889 value_t face_diameter = parent.compute_diameter(face_index, dim - 1);
891 const coefficient_t modulus = parent.modulus;
892 coefficient_t face_coefficient =
893 (k & 1 ? -1 + modulus : 1) * parent.get_coefficient(simplex) % modulus;
895 idx_below -= simplex_encoding(j, k + 1);
896 idx_above += simplex_encoding(j, k);
900 return parent.make_diameter_entry(face_diameter, face_index, face_coefficient);
905#if BOOST_VERSION >= 108100
906template <
class Key,
class T,
class H,
class E>
using Hash_map = boost::unordered_flat_map<Key, T, H, E>;
908template <
class Key,
class T,
class H,
class E>
using Hash_map = boost::unordered_map<Key, T, H, E>;
912 using size_t =
typename Filtration::size_t;
913 using coefficient_t =
typename Filtration::coefficient_t;
914 using coefficient_storage_t =
typename Filtration::coefficient_storage_t;
915 using dimension_t =
typename Filtration::dimension_t;
916 using value_t =
typename Filtration::value_t;
917 using vertex_t =
typename Filtration::vertex_t;
918 using simplex_t =
typename Filtration::simplex_t;
919 using diameter_simplex_t =
typename Filtration::diameter_simplex_t;
920 using entry_t =
typename Filtration::entry_t;
921 using diameter_entry_t =
typename Filtration::diameter_entry_t;
922 using Entry_hash =
typename Filtration::Entry_hash;
923 using Equal_index =
typename Filtration::Equal_index;
924 using Simplex_boundary_enumerator =
typename Filtration::Simplex_boundary_enumerator;
925 using Simplex_coboundary_enumerator =
typename Filtration::Simplex_coboundary_enumerator;
926 template<
class T>
using Greater_diameter_or_smaller_index =
typename Filtration::template Greater_diameter_or_smaller_index<T>;
928 typedef Compressed_sparse_matrix_<diameter_entry_t> Compressed_sparse_matrix;
929 typedef Hash_map<entry_t, size_t, Entry_hash, Equal_index> entry_hash_map;
933 const dimension_t dim_max;
934 const coefficient_t modulus;
935 const std::vector<coefficient_storage_t> multiplicative_inverse_;
936 mutable std::vector<diameter_entry_t> cofacet_entries;
937 mutable std::vector<vertex_t> vertices;
938 Simplex_boundary_enumerator facets;
940 Simplex_coboundary_enumerator cofacets1, cofacets2;
942 coefficient_t multiplicative_inverse(coefficient_t c)
const {
943 return multiplicative_inverse_[c];
947 : filt(
std::move(_filt)), n(filt.num_vertices()),
948 dim_max(
std::min<vertex_t>(_dim_max, n - 2)),
950 multiplicative_inverse_(multiplicative_inverse_vector<coefficient_storage_t>(_modulus)),
951 facets(filt), cofacets1(filt), cofacets2(filt)
956 std::optional<diameter_entry_t> get_zero_pivot_facet(
const diameter_entry_t simplex,
const dimension_t dim) {
957 facets.set_simplex(simplex, dim);
959 std::optional<diameter_entry_t> facet = facets.next();
961 if (get_diameter(*facet) == get_diameter(simplex))
return *facet;
966 std::optional<diameter_entry_t> get_zero_pivot_cofacet(
const diameter_entry_t simplex,
const dimension_t dim) {
967 cofacets1.set_simplex(simplex, dim);
969 std::optional<diameter_entry_t> cofacet = cofacets1.next_raw();
971 if (get_diameter(*cofacet) == get_diameter(simplex))
return *cofacet;
979 std::optional<diameter_entry_t> get_zero_apparent_facet(
const diameter_entry_t simplex,
const dimension_t dim) {
980 std::optional<diameter_entry_t> facet = get_zero_pivot_facet(simplex, dim);
981 if (!facet)
return std::nullopt;
982 std::optional<diameter_entry_t> cofacet = get_zero_pivot_cofacet(*facet, dim - 1);
983 if (!cofacet || filt.get_index(*cofacet) != filt.get_index(simplex))
return std::nullopt;
987 std::optional<diameter_entry_t> get_zero_apparent_cofacet(
const diameter_entry_t simplex,
const dimension_t dim) {
988 std::optional<diameter_entry_t> cofacet = get_zero_pivot_cofacet(simplex, dim);
989 if (!cofacet)
return std::nullopt;
990 std::optional<diameter_entry_t> facet = get_zero_pivot_facet(*cofacet, dim + 1);
991 if (!facet || filt.get_index(*facet) != filt.get_index(simplex))
return std::nullopt;
995 bool is_in_zero_apparent_pair(
const diameter_entry_t simplex,
const dimension_t dim) {
996 return get_zero_apparent_cofacet(simplex, dim) || get_zero_apparent_facet(simplex, dim);
999 void assemble_columns_to_reduce(std::vector<diameter_simplex_t>& simplices,
1000 std::vector<diameter_simplex_t>& columns_to_reduce,
1001 entry_hash_map& pivot_column_index, dimension_t dim) {
1003#ifdef GUDHI_INDICATE_PROGRESS
1004 std::cerr << clear_line <<
"assembling columns" << std::flush;
1005 std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step;
1008 columns_to_reduce.clear();
1009 std::vector<diameter_simplex_t> next_simplices;
1011 for (diameter_simplex_t& simplex : simplices) {
1012 cofacets2.set_simplex(filt.make_diameter_entry(simplex, 1), dim - 1);
1015 std::optional<diameter_entry_t> cofacet = cofacets2.next(
false);
1016 if (!cofacet)
break;
1017#ifdef GUDHI_INDICATE_PROGRESS
1018 if (std::chrono::steady_clock::now() > next) {
1019 std::cerr << clear_line <<
"assembling " << next_simplices.size()
1020 <<
" columns (processing " << std::distance(&simplices[0], &simplex)
1021 <<
"/" << simplices.size() <<
" simplices)" << std::flush;
1022 next = std::chrono::steady_clock::now() + time_step;
1025 if (dim < dim_max) next_simplices.push_back({get_diameter(*cofacet), filt.get_index(*cofacet)});
1027 if (!is_in_zero_apparent_pair(*cofacet, dim) &&
1028 (pivot_column_index.find(get_entry(*cofacet)) == pivot_column_index.end()))
1029 columns_to_reduce.push_back({get_diameter(*cofacet), filt.get_index(*cofacet)});
1033 if (dim < dim_max) simplices.swap(next_simplices);
1035#ifdef GUDHI_INDICATE_PROGRESS
1036 std::cerr << clear_line <<
"sorting " << columns_to_reduce.size() <<
" columns"
1040 std::sort(columns_to_reduce.begin(), columns_to_reduce.end(),
1041 Greater_diameter_or_smaller_index<diameter_simplex_t>(filt));
1042#ifdef GUDHI_INDICATE_PROGRESS
1043 std::cerr << clear_line << std::flush;
1047 template<
class OutPair>
1048 void compute_dim_0_pairs(std::vector<diameter_simplex_t>& edges,
1049 std::vector<diameter_simplex_t>& columns_to_reduce, OutPair& output_pair) {
1050 Union_find<vertex_t> dset(n);
1052 edges = filt.get_edges();
1053 std::sort(edges.rbegin(), edges.rend(),
1054 Greater_diameter_or_smaller_index<diameter_simplex_t>(filt));
1055 std::vector<vertex_t> vertices_of_edge(2);
1056 for (
auto e : edges) {
1058 filt.get_simplex_vertices(filt.get_index(e), 1, n, vertices_of_edge.rbegin());
1059 vertex_t u = dset.find(vertices_of_edge[0]), v = dset.find(vertices_of_edge[1]);
1062 if (get_diameter(e) != 0)
1063 output_pair(0, get_diameter(e));
1065 }
else if ((dim_max > 0) && !get_zero_apparent_cofacet(filt.make_diameter_entry(e, 1), 1))
1066 columns_to_reduce.push_back(e);
1068 if (dim_max > 0) std::reverse(columns_to_reduce.begin(), columns_to_reduce.end());
1070 for (vertex_t i = 0; i < n; ++i)
1071 if (dset.find(i) == i) output_pair(0, std::numeric_limits<value_t>::infinity());
1074 template <
typename Column> std::optional<diameter_entry_t> pop_pivot(Column& column) {
1075 while(!column.empty()) {
1076 diameter_entry_t pivot = column.top();
1079 if (column.empty() || filt.get_index(column.top()) != filt.get_index(pivot))
return pivot;
1080 coefficient_t sum = (filt.get_coefficient(pivot) + filt.get_coefficient(column.top())) % modulus;
1085 filt.set_coefficient(pivot, sum);
1088 return std::nullopt;
1091 template <
typename Column> std::optional<diameter_entry_t> get_pivot(Column& column) {
1093 std::optional<diameter_entry_t> result = pop_pivot(column);
1094 if (result) column.push(*result);
1099 template <
typename Column>
1100 std::optional<diameter_entry_t> init_coboundary_and_get_pivot(
const diameter_entry_t simplex,
1101 Column& working_coboundary,
const dimension_t dim,
1102 entry_hash_map& pivot_column_index) {
1103 bool check_for_emergent_pair =
true;
1104 cofacet_entries.clear();
1105 cofacets2.set_simplex(simplex, dim);
1107 std::optional<diameter_entry_t> cofacet = cofacets2.next();
1108 if (!cofacet)
break;
1109 cofacet_entries.push_back(*cofacet);
1110 if (check_for_emergent_pair && (get_diameter(simplex) == get_diameter(*cofacet))) {
1111 if ((pivot_column_index.find(get_entry(*cofacet)) == pivot_column_index.end()) &&
1112 !get_zero_apparent_facet(*cofacet, dim + 1))
1114 check_for_emergent_pair =
false;
1117 for (
auto cofacet : cofacet_entries) working_coboundary.push(cofacet);
1118 return get_pivot(working_coboundary);
1123 template <
typename Column>
1124 void add_simplex_coboundary(
const diameter_entry_t simplex,
const dimension_t dim,
1125 Column& working_reduction_column, Column& working_coboundary) {
1126 working_reduction_column.push(simplex);
1127 cofacets1.set_simplex(simplex, dim);
1129 std::optional<diameter_entry_t> cofacet = cofacets1.next();
1130 if (!cofacet)
break;
1131 working_coboundary.push(*cofacet);
1136 template <
typename Column>
1137 void add_coboundary(Compressed_sparse_matrix& reduction_matrix,
1138 const std::vector<diameter_simplex_t>& columns_to_reduce,
1139 const size_t index_column_to_add,
const coefficient_t factor,
1140 const dimension_t dim, Column& working_reduction_column,
1141 Column& working_coboundary) {
1142 diameter_entry_t column_to_add = filt.make_diameter_entry(columns_to_reduce[index_column_to_add], factor);
1143 add_simplex_coboundary(column_to_add, dim, working_reduction_column, working_coboundary);
1145 for (diameter_entry_t simplex : reduction_matrix.subrange(index_column_to_add)) {
1146 filt.set_coefficient(simplex, filt.get_coefficient(simplex) * factor % modulus);
1147 add_simplex_coboundary(simplex, dim, working_reduction_column, working_coboundary);
1151 template<
class OutPair>
1152 void compute_pairs(
const std::vector<diameter_simplex_t>& columns_to_reduce,
1153 entry_hash_map& pivot_column_index,
const dimension_t dim, OutPair& output_pair) {
1154 Compressed_sparse_matrix reduction_matrix;
1155 Greater_diameter_or_smaller_index<diameter_entry_t> cmp(filt);
1156 Heap<diameter_entry_t, std::vector<diameter_entry_t>,
1157 Greater_diameter_or_smaller_index<diameter_entry_t>>
1158 working_reduction_column(cmp), working_coboundary(cmp);
1160#ifdef GUDHI_INDICATE_PROGRESS
1161 std::chrono::steady_clock::time_point next = std::chrono::steady_clock::now() + time_step;
1163 for (
size_t index_column_to_reduce = 0; index_column_to_reduce < columns_to_reduce.size();
1164 ++index_column_to_reduce) {
1166 diameter_entry_t column_to_reduce = filt.make_diameter_entry(columns_to_reduce[index_column_to_reduce], 1);
1167 value_t diameter = get_diameter(column_to_reduce);
1169 reduction_matrix.append_column();
1171 working_reduction_column.clear(); working_coboundary.clear();
1173 std::optional<diameter_entry_t> pivot = init_coboundary_and_get_pivot(
1174 column_to_reduce, working_coboundary, dim, pivot_column_index);
1178#ifdef GUDHI_INDICATE_PROGRESS
1179 if (std::chrono::steady_clock::now() > next) {
1180 std::cerr << clear_line <<
"reducing column " << index_column_to_reduce + 1
1181 <<
"/" << columns_to_reduce.size() <<
" (diameter " << diameter <<
")"
1183 next = std::chrono::steady_clock::now() + time_step;
1187 auto pair = pivot_column_index.find(get_entry(*pivot));
1188 if (pair != pivot_column_index.end()) {
1189 entry_t other_pivot = pair->first;
1190 size_t index_column_to_add = pair->second;
1191 coefficient_t factor =
1192 modulus - filt.get_coefficient(*pivot) *
1193 multiplicative_inverse(filt.get_coefficient(other_pivot)) %
1197 add_coboundary(reduction_matrix, columns_to_reduce, index_column_to_add,
1198 factor, dim, working_reduction_column, working_coboundary);
1200 pivot = get_pivot(working_coboundary);
1201 }
else if (std::optional<diameter_entry_t> e = get_zero_apparent_facet(*pivot, dim + 1); e) {
1202 filt.set_coefficient(*e, modulus - filt.get_coefficient(*e));
1204 add_simplex_coboundary(*e, dim, working_reduction_column, working_coboundary);
1206 pivot = get_pivot(working_coboundary);
1208 value_t death = get_diameter(*pivot);
1209 output_pair(diameter, death);
1210 pivot_column_index.insert({get_entry(*pivot), index_column_to_reduce});
1214 std::optional<diameter_entry_t> e = pop_pivot(working_reduction_column);
1216 GUDHI_assert(filt.get_coefficient(*e) > 0);
1217 reduction_matrix.push_back(*e);
1222 output_pair(diameter, std::numeric_limits<value_t>::infinity());
1227#ifdef GUDHI_INDICATE_PROGRESS
1228 std::cerr << clear_line << std::flush;
1234 template<
class OutDim,
class OutPair>
1235 void compute_barcodes(OutDim&& output_dim, OutPair&& output_pair) {
1236 std::vector<diameter_simplex_t> simplices, columns_to_reduce;
1239 compute_dim_0_pairs(simplices, columns_to_reduce, output_pair);
1241 for (dimension_t dim = 1; dim <= dim_max; ++dim) {
1242 entry_hash_map pivot_column_index(0, filt, filt);
1243 pivot_column_index.reserve(columns_to_reduce.size());
1246 compute_pairs(columns_to_reduce, pivot_column_index, dim, output_pair);
1249 assemble_columns_to_reduce(simplices, columns_to_reduce, pivot_column_index,
1261 typedef std::size_t size_t;
1262 typedef float value_t;
1263 typedef int8_t dimension_t;
1264 typedef int vertex_t;
1265 typedef Gudhi::numbers::uint128_t simplex_t;
1267 typedef uint16_t coefficient_storage_t;
1268 typedef uint_least32_t coefficient_t;
1269 static const bool use_coefficients =
false;
1278template<
class Params,
class SimplexEncoding,
class DistanceMatrix,
class OutDim,
class OutPair>
1279void help2(DistanceMatrix&& dist,
int dim_max,
typename DistanceMatrix::value_t threshold,
unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) {
1280 typedef Rips_filtration<DistanceMatrix, SimplexEncoding, Params> Filt;
1282 Filt filt(std::move(dist), dim_max, threshold, modulus);
1283 Pcoh(std::move(filt), dim_max, modulus).compute_barcodes(output_dim, output_pair);
1285template<
bool use_coefficients_,
class simplex_t_,
class value_t_>
struct TParams {
1287 typedef std::size_t size_t;
1288 typedef value_t_ value_t;
1289 typedef int8_t dimension_t;
1290 typedef int vertex_t;
1291 typedef simplex_t_ simplex_t;
1292 typedef uint16_t coefficient_storage_t;
1293 typedef uint_least32_t coefficient_t;
1294 static const bool use_coefficients = use_coefficients_;
1296template<
class value_t_>
struct TParams2 {
1297 typedef int vertex_t;
1298 typedef value_t_ value_t;
1301template<
bool use_coefficients,
class DistanceMatrix,
class OutDim,
class OutPair>
1302void help1(DistanceMatrix&& dist,
int dim_max,
typename DistanceMatrix::value_t threshold,
unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) {
1303 auto n = dist.size();
1306 if (dim_max > n - 2) dim_max = n - 2;
1307 int bits_per_vertex = log2up(n);
1308 int bits_for_coeff = log2up(modulus - 1);
1309 int bitfield_size = bits_per_vertex * (dim_max + 2) + bits_for_coeff;
1310 if (bitfield_size <= 64) {
1311 typedef TParams<use_coefficients, uint64_t, typename DistanceMatrix::value_t> P;
1312 help2<P, Bitfield_encoding<P>>(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair);
1313 }
else if (bitfield_size <= 128) {
1314 typedef TParams<use_coefficients, Gudhi::numbers::uint128_t, typename DistanceMatrix::value_t> P;
1315 help2<P, Bitfield_encoding<P>>(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair);
1317 typedef TParams<use_coefficients, Gudhi::numbers::uint128_t, typename DistanceMatrix::value_t> P;
1318 help2<P, Cns_encoding<P>>(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair);
1323template<
class DistanceMatrix,
class OutDim,
class OutPair>
1324void ripser(DistanceMatrix dist,
int dim_max,
typename DistanceMatrix::value_t threshold,
unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) {
1326 help1<false>(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair);
1328 help1<true >(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair);
1331template<
class DMParams,
class OutDim,
class OutPair>
1332void ripser_auto(Sparse_distance_matrix<DMParams> dist,
int dim_max,
typename DMParams::value_t threshold,
unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) {
1333 ripser(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair);
1335template<
class DMParams, Compressed_matrix_layout Layout,
class OutDim,
class OutPair>
1336void ripser_auto(Compressed_distance_matrix<DMParams, Layout> dist,
int dim_max,
typename DMParams::value_t threshold,
unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) {
1337 typedef typename DMParams::value_t value_t;
1338 typedef typename DMParams::vertex_t vertex_t;
1339 if (threshold < std::numeric_limits<value_t>::max()) {
1340 Sparse_distance_matrix<DMParams> new_dist(dist, threshold);
1341 ripser(std::move(new_dist), dim_max, threshold, modulus, output_dim, output_pair);
1343 for (vertex_t i = 0; i < dist.size(); ++i) {
1344 value_t r_i = -std::numeric_limits<value_t>::infinity();
1345 for (vertex_t j = 0; j < dist.size(); ++j)
1346 r_i = std::max(r_i, dist(i, j));
1347 threshold = std::min(threshold, r_i);
1350 ripser(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair);
1354template<
class DistanceMatrix,
class OutDim,
class OutPair>
1355void ripser_auto(DistanceMatrix dist,
int dim_max,
typename DistanceMatrix::value_t threshold,
unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) {
1356 typedef typename DistanceMatrix::value_t value_t;
1357 typedef TParams2<value_t> P;
1358 if (threshold < std::numeric_limits<value_t>::max()) {
1359 Sparse_distance_matrix<P> new_dist(dist, threshold);
1360 ripser(std::move(new_dist), dim_max, threshold, modulus, output_dim, output_pair);
1362 Compressed_distance_matrix<P, LOWER_TRIANGULAR> new_dist(dist);
1363 ripser(std::move(new_dist), dim_max, threshold, modulus, output_dim, output_pair);
1367template<
class DistanceMatrix,
class OutDim,
class OutPair>
1368void ripser_auto(DistanceMatrix dist,
int dim_max,
typename DistanceMatrix::value_t threshold,
unsigned modulus, OutDim&& output_dim, OutPair&& output_pair) {
1369 typedef typename DistanceMatrix::vertex_t vertex_t;
1370 typedef typename DistanceMatrix::value_t value_t;
1371 typedef TParams2<value_t> P;
1372 if constexpr (std::is_same_v<typename DistanceMatrix::Category, Tag_sparse>) {
1373 ripser(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair);
1374 }
else if (threshold < std::numeric_limits<value_t>::max()) {
1375 Sparse_distance_matrix<P> new_dist(dist, threshold);
1376 ripser(std::move(new_dist), dim_max, threshold, modulus, output_dim, output_pair);
1377 }
else if constexpr (std::is_same_v<typename DistanceMatrix::Category, Tag_dense>) {
1378 for (vertex_t i = 0; i < dist.size(); ++i) {
1379 value_t r_i = -std::numeric_limits<value_t>::infinity();
1380 for (vertex_t j = 0; j < dist.size(); ++j)
1381 r_i = std::max(r_i, dist(i, j));
1382 threshold = std::min(threshold, r_i);
1385 ripser(std::move(dist), dim_max, threshold, modulus, output_dim, output_pair);
1387 Compressed_distance_matrix<P, LOWER_TRIANGULAR> new_dist(dist);
1388 ripser_auto(std::move(new_dist), dim_max, threshold, modulus, output_dim, output_pair);
Computes the persistent cohomology of a filtered complex.
Definition Persistent_cohomology.h:54