68#ifdef USE_ROBINHOOD_HASHMAP
70#include "robin_hood.h"
72template <
class Key,
class T,
class H,
class E>
73using hash_map = robin_hood::unordered_map<Key, T, H, E>;
75using hash = robin_hood::hash<Key>;
79template <
class Key,
class T,
class H,
class E>
80using hash_map = std::unordered_map<Key, T, H, E>;
82using hash = std::hash<Key>;
86#ifdef INDICATE_PROGRESS
87static const std::chrono::milliseconds time_step(40);
90static const std::string clear_line(
"\r\033[K");
92static const size_t num_coefficient_bits = 8;
95static const index_t max_simplex_index
96 = (__int128(1) << (8 *
sizeof(
index_t) - 1 - num_coefficient_bits)) - 1;
106#ifdef USE_COEFFICIENTS
107 (i > max_simplex_index)
111 throw std::overflow_error(
"simplex index " +
std::to_string((uint64_t)i)
112 +
" in filtration is larger than maximum index "
118 std::vector<index_t> B;
124 for(
index_t i = 0; i <= n; ++i) {
126 for(
index_t j = 1; j < std::min(i, k + 1); ++j)
128 = B[(i - 1) * offset + j - 1] + B[(i - 1) * offset + j];
130 B[i * offset + i] = 1;
136 assert(n <
index_t(B.size() / offset) && k <
index_t(offset) && n >= k - 1);
137 return B[n * offset + k];
145 return (modulus == 2) ? val & 1 : val % modulus;
149 return n > modulus / 2 ? n - modulus : n;
152std::vector<coefficient_t>
154 std::vector<coefficient_t> inverse(m);
160 inverse[a] = m - (inverse[m % a] * (m / a)) % m;
164#ifdef USE_COEFFICIENTS
168#define PACK(...) __pragma(pack(push, 1)) __VA_ARGS__ __pragma(pack(pop))
170#define PACK(...) __attribute__((__packed__)) __VA_ARGS__
175 index_t coefficient : num_coefficient_bits;
177 : index(_index), coefficient(_coefficient) {
181 entry_t() : index(0), coefficient(0) {
186 "size of entry_t is not the same as index_t");
201 return e.coefficient;
302template <
typename Entry>
316template <compressed_matrix_layout Layout>
329 template <
typename DistanceMatrix>
334 for(
size_t i = 1; i <
size(); ++i)
335 for(
size_t j = 0; j < i; ++j)
336 rows[i][j] = mat(i, j);
355 for(
size_t i = 1; i <
size(); ++i) {
364 for(
size_t i = 0; i <
size() - 1; ++i) {
366 pointer +=
size() - i - 2;
373 return i == j ? 0 : i < j ?
rows[j][i] :
rows[i][j];
379 return i == j ? 0 : i > j ?
rows[j][i] :
rows[i][j];
386 mutable std::vector<std::vector<index_diameter_t>::const_reverse_iterator>
388 mutable std::vector<std::vector<index_diameter_t>::const_reverse_iterator>
392 std::vector<std::vector<index_diameter_t>> &&_neighbors,
index_t _num_edges)
396 template <
typename DistanceMatrix>
399 for(
size_t i = 0; i <
size(); ++i)
400 for(
size_t j = 0; j <
size(); ++j)
401 if(i != j && mat(i, j) <= threshold) {
413 std::vector<std::vector<value_t>>
points;
422 return std::sqrt(std::inner_product(
424 std::plus<value_t>(),
434 std::vector<index_t> parent;
435 std::vector<uint8_t> rank;
445 while((z = parent[y]) != y)
447 while((z = parent[x]) != y) {
457 if(rank[x] > rank[y])
461 if(rank[x] == rank[y])
472T
end(std::pair<T, T> &p) {
476template <
typename ValueType>
478 std::vector<size_t> bounds;
479 std::vector<ValueType> entries;
481 using iterator =
typename std::vector<ValueType>::iterator;
482 using iterator_pair = std::pair<iterator, iterator>;
486 return bounds.size();
490 return {entries.begin() + (index == 0 ? 0 : bounds[index - 1]),
491 entries.begin() + bounds[index]};
495 bounds.push_back(entries.size());
500 entries.push_back(e);
511template <
typename DistanceMatrix>
513 const DistanceMatrix dist;
519 const std::vector<coefficient_t> multiplicative_inverse;
520 mutable std::vector<diameter_entry_t> cofacet_entries;
523 std::size_t operator()(
const entry_t &e)
const {
524#if defined(USE_ROBINHOOD_HASHMAP)
525 return robin_hood::hash<index_t>()(
::get_index(e));
546 : dist(
std::move(_dist)), n(dist.size()), dim_max(_dim_max),
547 threshold(_threshold), ratio(_ratio), modulus(_modulus),
548 binomial_coeff(n, dim_max + 2),
556 if(binomial_coeff(top, k) > idx) {
563 if(binomial_coeff(mid, k) > idx) {
574 return binomial_coeff(i, 2) + j;
577 template <
typename OutputIterator>
581 OutputIterator out)
const {
583 for(
index_t k = dim + 1; k > 0; --k) {
586 idx -= binomial_coeff(n_, k);
591 class simplex_coboundary_enumerator;
595 std::vector<diameter_index_t> &columns_to_reduce,
596 entry_hash_map &pivot_column_index,
598#ifdef INDICATE_PROGRESS
599 std::cerr << clear_line <<
"assembling columns" << std::flush;
600 std::chrono::steady_clock::time_point next
601 = std::chrono::steady_clock::now() + time_step;
604 columns_to_reduce.clear();
605 std::vector<diameter_index_t> next_simplices;
612#ifdef INDICATE_PROGRESS
613 if(std::chrono::steady_clock::now() > next) {
614 std::cerr << clear_line <<
"assembling " << next_simplices.size()
615 <<
" columns (processing "
616 << std::distance(&simplices[0], &simplex) <<
"/"
617 << simplices.size() <<
" simplices)" << std::flush;
618 next = std::chrono::steady_clock::now() + time_step;
621 auto cofacet = cofacets.
next();
624 next_simplices.emplace_back(
627 if(pivot_column_index.find(
get_entry(cofacet))
628 == pivot_column_index.end())
629 columns_to_reduce.emplace_back(
635 simplices.swap(next_simplices);
637#ifdef INDICATE_PROGRESS
638 std::cerr << clear_line <<
"sorting " << columns_to_reduce.size()
639 <<
" columns" << std::flush;
642 std::sort(columns_to_reduce.begin(), columns_to_reduce.end(),
645#ifdef INDICATE_PROGRESS
646 std::cerr << clear_line << std::flush;
651 std::vector<diameter_index_t> &columns_to_reduce,
652 std::vector<std::vector<pers_pair_t>> &ph) {
656 std::sort(edges.rbegin(), edges.rend(),
658 std::vector<index_t> vertices_of_edge(2);
659 for(
auto e : edges) {
662 v = dset.
find(vertices_of_edge[1]);
669 {u == dset.
find(vertices_of_edge[0]) ? vertices_of_edge[1]
670 : vertices_of_edge[0]},
673 {vertices_of_edge[0], vertices_of_edge[1]},
get_diameter(e)});
675 columns_to_reduce.push_back(e);
677 std::reverse(columns_to_reduce.begin(), columns_to_reduce.end());
679 for(
index_t i = 0; i < n; ++i) {
680 if(dset.
find(i) == i)
687 template <
typename Column>
690#ifdef USE_COEFFICIENTS
691 while(!column.empty()) {
693 pivot = column.top();
704 while(!column.empty()) {
705 pivot = column.top();
715 template <
typename Column>
723 template <
typename Column>
726 Column &working_coboundary,
728 entry_hash_map &pivot_column_index) {
729 bool check_for_emergent_pair =
true;
730 cofacet_entries.clear();
735 cofacet_entries.push_back(cofacet);
736 if(check_for_emergent_pair
738 if(pivot_column_index.find(
get_entry(cofacet))
739 == pivot_column_index.end())
741 check_for_emergent_pair =
false;
745 for(
auto cofacet : cofacet_entries)
746 working_coboundary.push(cofacet);
750 template <
typename Column>
753 Column &working_reduction_column,
754 Column &working_coboundary) {
755 working_reduction_column.push(simplex);
760 working_coboundary.push(cofacet);
764 template <
typename Column>
767 const std::vector<diameter_index_t> &columns_to_reduce,
768 const size_t index_column_to_add,
771 Column &working_reduction_column,
772 Column &working_coboundary) {
774 columns_to_reduce[index_column_to_add], factor);
776 column_to_add, dim, working_reduction_column, working_coboundary);
779 reduction_matrix.
subrange(index_column_to_add)) {
782 simplex, dim, working_reduction_column, working_coboundary);
788 std::vector<diameter_entry_t>,
792 entry_hash_map &pivot_column_index,
794 std::vector<std::vector<pers_pair_t>> &ph) {
796 size_t index_column_to_add;
798#ifdef INDICATE_PROGRESS
799 std::chrono::steady_clock::time_point next
800 = std::chrono::steady_clock::now() + time_step;
803 for(
size_t index_column_to_reduce = 0;
804 index_column_to_reduce < columns_to_reduce.size();
805 ++index_column_to_reduce) {
807 columns_to_reduce[index_column_to_reduce], 1);
815 working_reduction_column.push(column_to_reduce);
818 column_to_reduce, working_coboundary, dim, pivot_column_index);
821#ifdef INDICATE_PROGRESS
822 if(std::chrono::steady_clock::now() > next) {
823 std::cerr << clear_line <<
"reducing column "
824 << index_column_to_reduce + 1 <<
"/"
825 << columns_to_reduce.size() <<
" (diameter " << diameter
826 <<
")" << std::flush;
827 next = std::chrono::steady_clock::now() + time_step;
831 auto pair = pivot_column_index.find(
get_entry(pivot));
832 if(pair != pivot_column_index.end()) {
833 entry_t other_pivot = pair->first;
834 index_column_to_add = pair->second;
843 index_column_to_add, factor, dim,
844 working_reduction_column, working_coboundary);
849 if(death > diameter * ratio) {
850 std::vector<index_t> vertices_birth(dim + 1),
851 vertices_death(dim + 2);
853 get_index(column_to_reduce), dim, n, vertices_birth.rbegin());
855 get_index(pivot), dim + 1, n, vertices_death.rbegin());
860 pivot_column_index.insert(
861 {
get_entry(pivot), index_column_to_reduce});
875 std::vector<index_t> vertices_birth(dim + 1);
877 get_index(column_to_reduce), dim, n, vertices_birth.rbegin());
878 ph[dim].emplace_back(
885#ifdef INDICATE_PROGRESS
886 std::cerr << clear_line << std::flush;
892 std::vector<diameter_index_t> simplices, columns_to_reduce;
900 for(
index_t dim = 1; dim <= dim_max; ++dim) {
901 entry_hash_map pivot_column_index;
902 pivot_column_index.reserve(columns_to_reduce.size());
904 compute_pairs(columns_to_reduce, pivot_column_index, dim, ph);
908 simplices, columns_to_reduce, pivot_column_index, dim + 1);
916 index_t idx_below, idx_above, v, k;
917 std::vector<index_t> vertices;
928 : idx_below(
get_index(_simplex)), idx_above(0), v(parent.n - 1),
929 k(_dim + 1), vertices(_dim + 1), simplex(_simplex),
930 modulus(parent.modulus), dist(parent.dist),
931 binomial_coeff(parent.binomial_coeff) {
933 get_index(_simplex), _dim, parent.n, vertices.begin());
937 return (v >= k && (all_cofacets || binomial_coeff(v, k) > idx_below));
941 while((binomial_coeff(v, k) <= idx_below)) {
942 idx_below -= binomial_coeff(v, k);
943 idx_above += binomial_coeff(v, k + 1);
950 cofacet_diameter = std::max(cofacet_diameter, dist(v, w));
951 index_t cofacet_index = idx_above + binomial_coeff(v--, k + 1) + idx_below;
955 cofacet_diameter, cofacet_index, cofacet_coefficient);
961 index_t idx_below, idx_above, k;
962 std::vector<index_t> vertices;
967 std::vector<std::vector<index_diameter_t>::const_reverse_iterator>
969 std::vector<std::vector<index_diameter_t>::const_reverse_iterator>
977 : idx_below(
get_index(_simplex)), idx_above(0), k(_dim + 1),
978 vertices(_dim + 1), simplex(_simplex), modulus(parent.modulus),
979 dist(parent.dist), binomial_coeff(parent.binomial_coeff),
980 neighbor_it(dist.neighbor_it), neighbor_end(dist.neighbor_end) {
982 neighbor_end.clear();
986 for(
auto v : vertices) {
987 neighbor_it.push_back(dist.
neighbors[v].rbegin());
988 neighbor_end.push_back(dist.
neighbors[v].rend());
993 for(
auto &it0 = neighbor_it[0], &end0 = neighbor_end[0]; it0 != end0;
996 for(
size_t idx = 1; idx < neighbor_it.size(); ++idx) {
997 auto &it = neighbor_it[idx],
end = neighbor_end[idx];
1002 goto continue_outer;
1004 neighbor = std::max(neighbor, *it);
1006 while(k > 0 && vertices[k - 1] >
get_index(neighbor)) {
1009 idx_below -= binomial_coeff(vertices[k - 1], k);
1010 idx_above += binomial_coeff(vertices[k - 1], k + 1);
1024 = idx_above + binomial_coeff(
get_index(neighbor), k + 1) + idx_below;
1028 cofacet_diameter, cofacet_index, cofacet_coefficient);
1033std::vector<diameter_index_t>
1035 std::vector<diameter_index_t> edges;
1036 std::vector<index_t> vertices(2);
1037 for(
index_t index = binomial_coeff(n, 2); index-- > 0;) {
1038 get_simplex_vertices(index, 1, dist.size(), vertices.rbegin());
1039 value_t length = dist(vertices[0], vertices[1]);
1040 if(length <= threshold)
1041 edges.emplace_back(length, index);
1048 std::vector<diameter_index_t> edges;
1049 for(
index_t i = 0; i < n; ++i)
1050 for(
auto n_ : dist.neighbors[i]) {
1053 edges.emplace_back(
get_diameter(n_), get_edge_index(i, j));
1061 bool distanceMatrix,
1062 std::vector<std::vector<pers_pair_t>> &ph) {
1066 ph = std::vector<std::vector<pers_pair_t>>(
1067 dim_max + 1, std::vector<pers_pair_t>(0));
1069 if(!distanceMatrix) {
1073 std::move(dist), dim_max, threshold, ratio, modulus);
1074 ripser.compute_barcodes(ph);
1079 std::move(dist), dim_max, threshold, ratio, modulus);
1080 ripser.compute_barcodes(ph);
simplex_coboundary_enumerator(const diameter_entry_t _simplex, index_t _dim, const Ripser< compressed_lower_distance_matrix > &parent)
simplex_coboundary_enumerator(const diameter_entry_t _simplex, const index_t _dim, const Ripser< sparse_distance_matrix > &parent)
bool has_next(bool all_cofacets=true)
index_t get_edge_index(const index_t i, const index_t j) const
Ripser(DistanceMatrix &&_dist, index_t _dim_max, value_t _threshold, float _ratio, coefficient_t _modulus)
OutputIterator get_simplex_vertices(index_t idx, const index_t dim, index_t n_, OutputIterator out) const
void assemble_columns_to_reduce(std::vector< diameter_index_t > &simplices, std::vector< diameter_index_t > &columns_to_reduce, entry_hash_map &pivot_column_index, index_t dim)
diameter_entry_t init_coboundary_and_get_pivot(const diameter_entry_t simplex, Column &working_coboundary, const index_t &dim, entry_hash_map &pivot_column_index)
void compute_dim_0_pairs(std::vector< diameter_index_t > &edges, std::vector< diameter_index_t > &columns_to_reduce, std::vector< std::vector< pers_pair_t > > &ph)
void add_coboundary(compressed_sparse_matrix< diameter_entry_t > &reduction_matrix, const std::vector< diameter_index_t > &columns_to_reduce, const size_t index_column_to_add, const coefficient_t factor, const size_t &dim, Column &working_reduction_column, Column &working_coboundary)
void compute_barcodes(std::vector< std::vector< pers_pair_t > > &ph)
diameter_entry_t get_pivot(Column &column)
void add_simplex_coboundary(const diameter_entry_t simplex, const index_t &dim, Column &working_reduction_column, Column &working_coboundary)
std::vector< diameter_index_t > get_edges()
index_t get_max_vertex(const index_t idx, const index_t k, const index_t n_) const
void compute_pairs(std::vector< diameter_index_t > &columns_to_reduce, entry_hash_map &pivot_column_index, index_t dim, std::vector< std::vector< pers_pair_t > > &ph)
diameter_entry_t pop_pivot(Column &column)
std::priority_queue< diameter_entry_t, std::vector< diameter_entry_t >, greater_diameter_or_smaller_index< diameter_entry_t > > working_t
binomial_coeff_table(index_t n, index_t k)
index_t operator()(index_t n, index_t k) const
std::vector< value_t > distances
value_t operator()(const index_t i, const index_t j) const
compressed_distance_matrix(std::vector< value_t > &&_distances)
std::vector< value_t * > rows
compressed_distance_matrix(const DistanceMatrix &mat)
void push_back(const ValueType e)
iterator_pair subrange(const index_t index)
union_find(const index_t n)
void link(index_t x, index_t y)
std::pair< simplex_t, value_t > simplex_diam_t
void ripser(std::vector< std::vector< value_t > > points, value_t threshold, index_t dim_max, bool distanceMatrix, std::vector< std::vector< pers_pair_t > > &ph)
std::string to_string(__int128)
std::ostream & operator<<(std::ostream &o, Node const &n)
std::pair< value_t, index_t > diameter_index_t
index_t get_coefficient(const entry_t &)
index_t get_index(const entry_t &i)
coefficient_t get_modulo(const coefficient_t val, const coefficient_t modulus)
T end(std::pair< T, T > &p)
value_t get_diameter(const diameter_index_t &i)
coefficient_t normalize(const coefficient_t n, const coefficient_t modulus)
std::pair< index_t, value_t > index_diameter_t
std::vector< coefficient_t > multiplicative_inverse_vector(const coefficient_t m)
entry_t make_entry(index_t _index, coefficient_t)
const entry_t & get_entry(const entry_t &e)
void set_coefficient(entry_t &, const coefficient_t)
void check_overflow(index_t i)
std::unordered_map< Key, T, H, E > hash_map
T begin(std::pair< T, T > &p)
diameter_entry_t(value_t _diameter, index_t _index, coefficient_t _coefficient)
diameter_entry_t()=default
diameter_entry_t(const index_t &_index)
diameter_entry_t(const diameter_index_t &_diameter_index, coefficient_t _coefficient)
value_t operator()(const index_t i, const index_t j) const
euclidean_distance_matrix(std::vector< std::vector< value_t > > &&_points)
std::vector< std::vector< value_t > > points
bool operator()(const Entry &a, const Entry &b)
std::vector< std::vector< index_diameter_t > > neighbors
std::vector< std::vector< index_diameter_t >::const_reverse_iterator > neighbor_end
sparse_distance_matrix(std::vector< std::vector< index_diameter_t > > &&_neighbors, index_t _num_edges)
std::vector< std::vector< index_diameter_t >::const_reverse_iterator > neighbor_it
sparse_distance_matrix(const DistanceMatrix &mat, const value_t threshold)