From 1aa3492247cf15590a2b4f4dd6098892b7520630 Mon Sep 17 00:00:00 2001 From: julian Date: Mon, 1 Mar 2021 15:06:49 +0100 Subject: [PATCH 1/5] [C++] Repalce USE_COEFFICIENTS with templates Signed-off-by: julian --- ripser/ripser.cpp | 385 ++++++++++++++++++++++++++++------------------ 1 file changed, 235 insertions(+), 150 deletions(-) diff --git a/ripser/ripser.cpp b/ripser/ripser.cpp index f2b7531..7620514 100644 --- a/ripser/ripser.cpp +++ b/ripser/ripser.cpp @@ -85,12 +85,7 @@ static const index_t max_simplex_index = void check_overflow(index_t i) { - if -#ifdef USE_COEFFICIENTS - (i > max_simplex_index) -#else - (i < 0) -#endif + if (i > max_simplex_index) throw std::overflow_error( "simplex index " + std::to_string((uint64_t) i) + " in filtration is larger than maximum index " + @@ -151,54 +146,63 @@ std::vector multiplicative_inverse_vector(const coefficient_t m) return inverse; } -#ifdef USE_COEFFICIENTS - // https://stackoverflow.com/a/3312896/13339777 #ifdef _MSC_VER -#define PACK( ... ) __pragma( pack(push, 1) ) __VA_ARGS__ __pragma( pack(pop)) +#define PACK(...) __pragma(pack(push, 1)) __VA_ARGS__ __pragma(pack(pop)) #else -#define PACK( ... ) __attribute__((__packed__)) __VA_ARGS__ +#define PACK(...) __attribute__((__packed__)) __VA_ARGS__ #endif -PACK(struct entry_t { +PACK(struct entry_t_N_ { index_t index : 8 * sizeof(index_t) - num_coefficient_bits; index_t coefficient : num_coefficient_bits; - entry_t(index_t _index, coefficient_t _coefficient) + entry_t_N_(index_t _index, coefficient_t _coefficient) : index(_index), coefficient(_coefficient) { } - entry_t(index_t _index) : index(_index), coefficient(0) {} - entry_t() : index(0), coefficient(0) {} + entry_t_N_(index_t _index) : index(_index), coefficient(0) {} + entry_t_N_() : index(0), coefficient(0) {} }); -static_assert(sizeof(entry_t) == sizeof(index_t), - "size of entry_t is not the same as index_t"); +using entry_t_2 = index_t; +using entry_t_N = struct entry_t_N_; +static_assert(sizeof(entry_t_N) == sizeof(index_t), + "size of entry_t_N is not the same as index_t"); -entry_t make_entry(index_t i, coefficient_t c) { return entry_t(i, c); } -index_t get_index(const entry_t& e) { return e.index; } -index_t get_coefficient(const entry_t& e) { return e.coefficient; } -void set_coefficient(entry_t& e, const coefficient_t c) { e.coefficient = c; } +const index_t get_index(const entry_t_2& i) { return i; } +const index_t get_index(const entry_t_N& i) { return i.index; } -std::ostream& operator<<(std::ostream& stream, const entry_t& e) +index_t get_coefficient(const entry_t_2& i) { return 1; } +index_t get_coefficient(const entry_t_N& i) { return i.coefficient; } + +template +typename std::enable_if::value, entry_t_2>::type +make_entry(index_t i, coefficient_t c) { - stream << get_index(e) << ":" << get_coefficient(e); - return stream; + return entry_t_2(i); } -#else - -typedef index_t entry_t; -const index_t get_index(const entry_t& i) { return i; } -index_t get_coefficient(const entry_t& i) { return 1; } -entry_t make_entry(index_t _index, coefficient_t _value) +template +typename std::enable_if::value, entry_t_N>::type +make_entry(index_t i, coefficient_t c) { - return entry_t(_index); + return entry_t_N(i, c); } -void set_coefficient(entry_t& e, const coefficient_t c) {} -#endif +void set_coefficient(entry_t_2& e, const coefficient_t c) {} +void set_coefficient(entry_t_N& e, const coefficient_t c) { e.coefficient = c; } -const entry_t& get_entry(const entry_t& e) { return e; } +std::ostream& operator<<(std::ostream& stream, const entry_t_N& e) +{ + stream << get_index(e) << ":" << get_coefficient(e); + return stream; +} + +template +const entry_t& get_entry(const entry_t& e) +{ + return e; +} typedef std::pair diameter_index_t; value_t get_diameter(const diameter_index_t& i) { return i.first; } @@ -208,35 +212,52 @@ typedef std::pair index_diameter_t; index_t get_index(const index_diameter_t& i) { return i.first; } value_t get_diameter(const index_diameter_t& i) { return i.second; } +template struct diameter_entry_t : std::pair { using std::pair::pair; diameter_entry_t() {} diameter_entry_t(value_t _diameter, index_t _index, coefficient_t _coefficient) - : diameter_entry_t(_diameter, make_entry(_index, _coefficient)) + : diameter_entry_t(_diameter, make_entry(_index, _coefficient)) { } diameter_entry_t(const diameter_index_t& _diameter_index, coefficient_t _coefficient) - : diameter_entry_t(get_diameter(_diameter_index), - make_entry(get_index(_diameter_index), _coefficient)) + : diameter_entry_t( + get_diameter(_diameter_index), + make_entry(get_index(_diameter_index), _coefficient)) { } diameter_entry_t(const index_t& _index) : diameter_entry_t(0, _index, 0) {} }; -const entry_t& get_entry(const diameter_entry_t& p) { return p.second; } -entry_t& get_entry(diameter_entry_t& p) { return p.second; } -const index_t get_index(const diameter_entry_t& p) +template +const entry_t& get_entry(const diameter_entry_t& p) +{ + return p.second; +} +template +entry_t& get_entry(diameter_entry_t& p) +{ + return p.second; +} +template +const index_t get_index(const diameter_entry_t& p) { return get_index(get_entry(p)); } -const coefficient_t get_coefficient(const diameter_entry_t& p) +template +const coefficient_t get_coefficient(const diameter_entry_t& p) { return get_coefficient(get_entry(p)); } -const value_t& get_diameter(const diameter_entry_t& p) { return p.first; } -void set_coefficient(diameter_entry_t& p, const coefficient_t c) +template +const value_t& get_diameter(const diameter_entry_t& p) +{ + return p.first; +} +template +void set_coefficient(diameter_entry_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); } @@ -508,7 +529,25 @@ typedef struct { int num_edges; } ripserResults; -template +template +class ripser; +template +class simplex_coboundary_enumerator; +template +std::vector +get_edges(const ripser& _parent); +template +std::vector +get_edges(const ripser& _parent); +template +value_t +get_vertex_birth(const ripser&, + index_t i); +template +value_t get_vertex_birth(const ripser&, + index_t i); + +template class ripser { const DistanceMatrix dist; @@ -518,7 +557,7 @@ class ripser const coefficient_t modulus; const binomial_coeff_table binomial_coeff; const std::vector multiplicative_inverse; - mutable std::vector cofacet_entries; + mutable std::vector> cofacet_entries; // If this flag is off, don't extract the representative cocycles to save // time const int do_cocycles; @@ -543,6 +582,23 @@ class ripser typedef hash_map entry_hash_map; + /* Make friend and class `friend` in order to access private data */ + template + friend class simplex_coboundary_enumerator; + template + friend std::vector + get_edges(const ripser& _parent); + template + friend std::vector + get_edges(const ripser& _parent); + template + friend value_t + get_vertex_birth(const ripser&, + index_t i); + template + friend value_t get_vertex_birth(const ripser&, + index_t i); + public: mutable std::vector> births_and_deaths_by_dim; mutable std::vector>> cocycles_by_dim; @@ -603,8 +659,6 @@ class ripser return out; } - class simplex_coboundary_enumerator; - void assemble_columns_to_reduce(std::vector& simplices, std::vector& columns_to_reduce, @@ -620,8 +674,8 @@ class ripser std::vector next_simplices; for (diameter_index_t& simplex : simplices) { - simplex_coboundary_enumerator cofacets(diameter_entry_t(simplex, 1), - dim, *this); + simplex_coboundary_enumerator cofacets( + diameter_entry_t(simplex, 1), dim, *this); while (cofacets.has_next(false)) { #ifdef INDICATE_PROGRESS @@ -662,7 +716,6 @@ class ripser #endif } - value_t get_vertex_birth(index_t i); void compute_dim_0_pairs(std::vector& edges, std::vector& columns_to_reduce) { @@ -670,10 +723,10 @@ class ripser // lower star) union_find dset(n); for (index_t i = 0; i < n; i++) { - dset.set_birth(i, get_vertex_birth(i)); + dset.set_birth(i, get_vertex_birth(*this, i)); } - edges = get_edges(); + edges = get_edges(*this); std::sort(edges.rbegin(), edges.rend(), greater_diameter_or_smaller_index()); std::vector vertices_of_edge(2); @@ -709,54 +762,56 @@ class ripser } template - diameter_entry_t pop_pivot(Column& column) + diameter_entry_t pop_pivot(Column& column) { - diameter_entry_t pivot(-1); -#ifdef USE_COEFFICIENTS - while (!column.empty()) { - if (get_coefficient(pivot) == 0) + diameter_entry_t pivot(-1); + if (std::is_same::value) { + while (!column.empty()) { + if (get_coefficient(pivot) == 0) + pivot = column.top(); + else if (get_index(column.top()) != get_index(pivot)) + return pivot; + else + set_coefficient(pivot, + get_modulo((get_coefficient(pivot) + + get_coefficient(column.top())), + modulus)); + column.pop(); + } + return (get_coefficient(pivot) == 0) ? -1 : pivot; + } else { + while (!column.empty()) { pivot = column.top(); - else if (get_index(column.top()) != get_index(pivot)) - return pivot; - else - set_coefficient(pivot, - get_modulo((get_coefficient(pivot) + - get_coefficient(column.top())), - modulus)); - column.pop(); - } - return (get_coefficient(pivot) == 0) ? -1 : pivot; -#else - while (!column.empty()) { - pivot = column.top(); - column.pop(); - if (column.empty() || get_index(column.top()) != get_index(pivot)) - return pivot; - column.pop(); + column.pop(); + if (column.empty() || + get_index(column.top()) != get_index(pivot)) + return pivot; + column.pop(); + } } return -1; -#endif } template - diameter_entry_t get_pivot(Column& column) + diameter_entry_t get_pivot(Column& column) { - diameter_entry_t result = pop_pivot(column); + diameter_entry_t result = pop_pivot(column); if (get_index(result) != -1) column.push(result); return result; } template - diameter_entry_t init_coboundary_and_get_pivot( - const diameter_entry_t simplex, Column& working_coboundary, + 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) { bool check_for_emergent_pair = true; cofacet_entries.clear(); - simplex_coboundary_enumerator cofacets(simplex, dim, *this); + simplex_coboundary_enumerator cofacets( + simplex, dim, *this); while (cofacets.has_next()) { - diameter_entry_t cofacet = cofacets.next(); + diameter_entry_t cofacet = cofacets.next(); if (get_diameter(cofacet) <= threshold) { cofacet_entries.push_back(cofacet); if (check_for_emergent_pair && @@ -774,34 +829,35 @@ class ripser } template - void add_simplex_coboundary(const diameter_entry_t simplex, + void add_simplex_coboundary(const diameter_entry_t simplex, const index_t& dim, Column& working_reduction_column, Column& working_coboundary) { working_reduction_column.push(simplex); - simplex_coboundary_enumerator cofacets(simplex, dim, *this); + simplex_coboundary_enumerator cofacets( + simplex, dim, *this); while (cofacets.has_next()) { - diameter_entry_t cofacet = cofacets.next(); + diameter_entry_t cofacet = cofacets.next(); if (get_diameter(cofacet) <= threshold) working_coboundary.push(cofacet); } } template - void - add_coboundary(compressed_sparse_matrix& reduction_matrix, - const std::vector& 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 add_coboundary( + compressed_sparse_matrix>& reduction_matrix, + const std::vector& 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) { - diameter_entry_t column_to_add(columns_to_reduce[index_column_to_add], - factor); + diameter_entry_t column_to_add( + columns_to_reduce[index_column_to_add], factor); add_simplex_coboundary(column_to_add, dim, working_reduction_column, working_coboundary); - for (diameter_entry_t simplex : + for (diameter_entry_t simplex : reduction_matrix.subrange(index_column_to_add)) { set_coefficient(simplex, get_coefficient(simplex) * factor % modulus); @@ -811,10 +867,10 @@ class ripser } using working_t = std::priority_queue< - diameter_entry_t, std::vector, - greater_diameter_or_smaller_index>; + diameter_entry_t, std::vector>, + greater_diameter_or_smaller_index>>; - diameter_entry_t cocycle_e; + diameter_entry_t cocycle_e; std::vector cocycle_simplex; std::vector thiscocycle; inline void compute_cocycles(working_t cocycle, index_t dim) @@ -837,7 +893,7 @@ class ripser void compute_pairs(std::vector& columns_to_reduce, entry_hash_map& pivot_column_index, index_t dim) { - compressed_sparse_matrix reduction_matrix; + compressed_sparse_matrix> reduction_matrix; size_t index_column_to_add; #ifdef INDICATE_PROGRESS @@ -848,7 +904,7 @@ class ripser for (size_t index_column_to_reduce = 0; index_column_to_reduce < columns_to_reduce.size(); ++index_column_to_reduce) { - diameter_entry_t column_to_reduce( + diameter_entry_t column_to_reduce( columns_to_reduce[index_column_to_reduce], 1); value_t diameter = get_diameter(column_to_reduce); @@ -859,7 +915,7 @@ class ripser working_reduction_column.push(column_to_reduce); - diameter_entry_t pivot = init_coboundary_and_get_pivot( + diameter_entry_t pivot = init_coboundary_and_get_pivot( column_to_reduce, working_coboundary, dim, pivot_column_index); while (true) { @@ -907,7 +963,7 @@ class ripser pop_pivot(working_reduction_column); while (true) { - diameter_entry_t e = + diameter_entry_t e = pop_pivot(working_reduction_column); if (get_index(e) == -1) @@ -935,7 +991,6 @@ class ripser #endif } - std::vector get_edges(); void compute_barcodes() { std::vector simplices, columns_to_reduce; @@ -962,35 +1017,38 @@ class ripser } }; -template <> -value_t ripser::get_vertex_birth(index_t i) +template +value_t +get_vertex_birth(const ripser&, + index_t) { // TODO: Dummy for now; nonzero vertex births are only done through // sparse matrices at the moment return 0.0; } -template <> -value_t ripser::get_vertex_birth(index_t i) +template +value_t get_vertex_birth(const ripser& parent_, + index_t i) { - return dist.vertex_births[i]; + return parent_.dist.vertex_births[i]; } -template <> -class ripser::simplex_coboundary_enumerator +template +class simplex_coboundary_enumerator { private: index_t idx_below, idx_above, v, k; std::vector vertices; - const diameter_entry_t simplex; + const diameter_entry_t simplex; const coefficient_t modulus; const compressed_lower_distance_matrix& dist; const binomial_coeff_table& binomial_coeff; public: simplex_coboundary_enumerator( - const diameter_entry_t _simplex, index_t _dim, - const ripser& parent) + const diameter_entry_t _simplex, index_t _dim, + const ripser& parent) : idx_below(get_index(_simplex)), idx_above(0), v(parent.n - 1), k(_dim + 1), vertices(_dim + 1), simplex(_simplex), modulus(parent.modulus), dist(parent.dist), @@ -1005,7 +1063,7 @@ class ripser::simplex_coboundary_enumerator return (v >= k && (all_cofacets || binomial_coeff(v, k) > idx_below)); } - diameter_entry_t next() + diameter_entry_t next() { while ((binomial_coeff(v, k) <= idx_below)) { idx_below -= binomial_coeff(v, k); @@ -1021,17 +1079,17 @@ class ripser::simplex_coboundary_enumerator idx_above + binomial_coeff(v--, k + 1) + idx_below; coefficient_t cofacet_coefficient = (k & 1 ? modulus - 1 : 1) * get_coefficient(simplex) % modulus; - return diameter_entry_t(cofacet_diameter, cofacet_index, - cofacet_coefficient); + return diameter_entry_t(cofacet_diameter, cofacet_index, + cofacet_coefficient); } }; -template <> -class ripser::simplex_coboundary_enumerator +template +class simplex_coboundary_enumerator { index_t idx_below, idx_above, k; std::vector vertices; - const diameter_entry_t simplex; + const diameter_entry_t simplex; const coefficient_t modulus; const sparse_distance_matrix& dist; const binomial_coeff_table& binomial_coeff; @@ -1042,9 +1100,9 @@ class ripser::simplex_coboundary_enumerator index_diameter_t neighbor; public: - simplex_coboundary_enumerator(const diameter_entry_t _simplex, - const index_t _dim, - const ripser& parent) + simplex_coboundary_enumerator( + const diameter_entry_t _simplex, const index_t _dim, + const ripser& parent) : idx_below(get_index(_simplex)), idx_above(0), k(_dim + 1), vertices(_dim + 1), simplex(_simplex), modulus(parent.modulus), dist(parent.dist), binomial_coeff(parent.binomial_coeff), @@ -1090,7 +1148,7 @@ class ripser::simplex_coboundary_enumerator return false; } - diameter_entry_t next() + diameter_entry_t next() { ++neighbor_it[0]; value_t cofacet_diameter = @@ -1099,35 +1157,38 @@ class ripser::simplex_coboundary_enumerator idx_above + binomial_coeff(get_index(neighbor), k + 1) + idx_below; coefficient_t cofacet_coefficient = (k & 1 ? modulus - 1 : 1) * get_coefficient(simplex) % modulus; - return diameter_entry_t(cofacet_diameter, cofacet_index, - cofacet_coefficient); + return diameter_entry_t(cofacet_diameter, cofacet_index, + cofacet_coefficient); } }; -template <> +template std::vector -ripser::get_edges() +get_edges(const ripser& _parent) { std::vector edges; std::vector vertices(2); - for (index_t index = binomial_coeff(n, 2); index-- > 0;) { - get_simplex_vertices(index, 1, dist.size(), vertices.rbegin()); - value_t length = dist(vertices[0], vertices[1]); - if (length <= threshold) + for (index_t index = _parent.binomial_coeff(_parent.n, 2); index-- > 0;) { + _parent.get_simplex_vertices(index, 1, _parent.dist.size(), + vertices.rbegin()); + value_t length = _parent.dist(vertices[0], vertices[1]); + if (length <= _parent.threshold) edges.push_back({length, index}); } return edges; } -template <> -std::vector ripser::get_edges() +template +std::vector +get_edges(const ripser& _parent) { std::vector edges; - for (index_t i = 0; i < n; ++i) - for (auto n : dist.neighbors[i]) { + for (index_t i = 0; i < _parent.n; ++i) + for (auto n : _parent.dist.neighbors[i]) { index_t j = get_index(n); if (i > j) - edges.push_back({get_diameter(n), get_edge_index(i, j)}); + edges.push_back( + {get_diameter(n), _parent.get_edge_index(i, j)}); } return edges; } @@ -1174,16 +1235,31 @@ ripserResults rips_dm(float* D, int N, int modulus, int dim_max, ripserResults res; if (threshold >= max) { - ripser r( + if (modulus == 2) { + ripser r( + std::move(dist), dim_max, threshold, ratio, modulus, do_cocycles); + r.compute_barcodes(); + r.copy_results(res); + } else { + ripser r( std::move(dist), dim_max, threshold, ratio, modulus, do_cocycles); r.compute_barcodes(); r.copy_results(res); + } } else { - ripser r( + if (modulus == 2) { + ripser r( sparse_distance_matrix(std::move(dist), threshold), dim_max, threshold, ratio, modulus, do_cocycles); r.compute_barcodes(); r.copy_results(res); + } else { + ripser r( + sparse_distance_matrix(std::move(dist), threshold), dim_max, + threshold, ratio, modulus, do_cocycles); + r.compute_barcodes(); + r.copy_results(res); + } } res.num_edges = num_edges; return res; @@ -1195,11 +1271,21 @@ ripserResults rips_dm_sparse(int* I, int* J, float* V, int NEdges, int N, { // TODO: This seems like a dummy parameter at the moment float ratio = 1.0; + ripserResults res; // Setup distance matrix and figure out threshold - ripser r( - sparse_distance_matrix(I, J, V, NEdges, N, threshold), dim_max, - threshold, ratio, modulus, do_cocycles); - r.compute_barcodes(); + if (modulus == 2) { + ripser r( + sparse_distance_matrix(I, J, V, NEdges, N, threshold), dim_max, + threshold, ratio, modulus, do_cocycles); + r.compute_barcodes(); + r.copy_results(res); + } else { + ripser r( + sparse_distance_matrix(I, J, V, NEdges, N, threshold), dim_max, + threshold, ratio, modulus, do_cocycles); + r.compute_barcodes(); + r.copy_results(res); + } // Report the number of edges that were added int num_edges = 0; for (int idx = 0; idx < NEdges; idx++) { @@ -1207,8 +1293,6 @@ ripserResults rips_dm_sparse(int* I, int* J, float* V, int NEdges, int N, num_edges++; } } - ripserResults res; - r.copy_results(res); res.num_edges = num_edges; return res; } @@ -1274,24 +1358,25 @@ extern "C" { /* C interface to Ripser. - Results are passed through output arguments. The arrays are allocated in this - function and have to be freed manually by the caller. + Results are passed through output arguments. The arrays are allocated in + this function and have to be freed manually by the caller. Output arguments: * n_intervals: number of intervals per dimension. (length = dim_max + 1) * births_and_deaths: births and deaths of all dimension in a flat array. (length = 2 * sum(n_intervals)) - * cocycle_length: lengths of individual cocycles. (length = sum(n_intervals)) - * cocycles: cocycles stored in a flat array. (length = sum(cocycle_length)) - Input arguments: + * cocycle_length: lengths of individual cocycles. (length = + sum(n_intervals)) + * cocycles: cocycles stored in a flat array. (length = + sum(cocycle_length)) Input arguments: * D: lower triangle of the distance matrix in a flat array. * N: length of D. - * modulus: Compute homology with coefficients in the prime field Z/pZ. p must - be a prime number. + * modulus: Compute homology with coefficients in the prime field Z/pZ. p + must be a prime number. * dim_max: Compute persistent homology up to this dimension * threshold: Compute Rips complexes up to this diameter - * do_cocycles: If nonzero, calculate cocycles and write them to cocycle_length - and cocycles. + * do_cocycles: If nonzero, calculate cocycles and write them to + cocycle_length and cocycles. Returns number of edges. */ From ba3c335f2ae4dec8c8448853077298fe8bc16a22 Mon Sep 17 00:00:00 2001 From: julian Date: Mon, 1 Mar 2021 15:48:17 +0100 Subject: [PATCH 2/5] [C++] Refactor call of ripser with different template values Signed-off-by: julian --- ripser/ripser.cpp | 76 +++++++++++++++++++++-------------------------- 1 file changed, 34 insertions(+), 42 deletions(-) diff --git a/ripser/ripser.cpp b/ripser/ripser.cpp index 7620514..e359a47 100644 --- a/ripser/ripser.cpp +++ b/ripser/ripser.cpp @@ -1193,6 +1193,32 @@ get_edges(const ripser& _parent) return edges; } +template +ripserResults get_ripser_modulus(DistanceMatrix&& dist, int modulus, + int dim_max, float threshold, int do_cocycles) +{ + // TODO: This seems like a dummy parameter at the moment + float ratio = 1.0; + ripserResults res; + + switch (modulus) { + case 2: { + ripser r(std::move(dist), dim_max, threshold, + ratio, modulus, do_cocycles); + r.compute_barcodes(); + r.copy_results(res); + break; + } + default: { + ripser r(std::move(dist), dim_max, threshold, + ratio, modulus, do_cocycles); + r.compute_barcodes(); + r.copy_results(res); + } + } + return res; +} + ripserResults rips_dm(float* D, int N, int modulus, int dim_max, float threshold, int do_cocycles) { @@ -1201,9 +1227,6 @@ ripserResults rips_dm(float* D, int N, int modulus, int dim_max, compressed_lower_distance_matrix dist = compressed_lower_distance_matrix( compressed_upper_distance_matrix(std::move(distances))); - // TODO: This seems like a dummy parameter at the moment - float ratio = 1.0; - value_t min = std::numeric_limits::infinity(), max = -std::numeric_limits::infinity(), max_finite = max; int num_edges = 0; @@ -1235,31 +1258,12 @@ ripserResults rips_dm(float* D, int N, int modulus, int dim_max, ripserResults res; if (threshold >= max) { - if (modulus == 2) { - ripser r( - std::move(dist), dim_max, threshold, ratio, modulus, do_cocycles); - r.compute_barcodes(); - r.copy_results(res); - } else { - ripser r( - std::move(dist), dim_max, threshold, ratio, modulus, do_cocycles); - r.compute_barcodes(); - r.copy_results(res); - } + res = get_ripser_modulus(std::move(dist), modulus, dim_max, threshold, + do_cocycles); } else { - if (modulus == 2) { - ripser r( - sparse_distance_matrix(std::move(dist), threshold), dim_max, - threshold, ratio, modulus, do_cocycles); - r.compute_barcodes(); - r.copy_results(res); - } else { - ripser r( - sparse_distance_matrix(std::move(dist), threshold), dim_max, - threshold, ratio, modulus, do_cocycles); - r.compute_barcodes(); - r.copy_results(res); - } + res = get_ripser_modulus( + sparse_distance_matrix(std::move(dist), threshold), modulus, + dim_max, threshold, do_cocycles); } res.num_edges = num_edges; return res; @@ -1269,23 +1273,11 @@ ripserResults rips_dm_sparse(int* I, int* J, float* V, int NEdges, int N, int modulus, int dim_max, float threshold, int do_cocycles) { - // TODO: This seems like a dummy parameter at the moment - float ratio = 1.0; ripserResults res; // Setup distance matrix and figure out threshold - if (modulus == 2) { - ripser r( - sparse_distance_matrix(I, J, V, NEdges, N, threshold), dim_max, - threshold, ratio, modulus, do_cocycles); - r.compute_barcodes(); - r.copy_results(res); - } else { - ripser r( - sparse_distance_matrix(I, J, V, NEdges, N, threshold), dim_max, - threshold, ratio, modulus, do_cocycles); - r.compute_barcodes(); - r.copy_results(res); - } + res = get_ripser_modulus( + sparse_distance_matrix(I, J, V, NEdges, N, threshold), modulus, dim_max, + threshold, do_cocycles); // Report the number of edges that were added int num_edges = 0; for (int idx = 0; idx < NEdges; idx++) { From f0a4fd34f75d6a19b48d4c1e941098d92e6ba6d6 Mon Sep 17 00:00:00 2001 From: julian Date: Mon, 1 Mar 2021 15:56:29 +0100 Subject: [PATCH 3/5] Remove unused defines from setup.py Signed-off-by: julian --- setup.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/setup.py b/setup.py index 5f44f2f..41828c0 100644 --- a/setup.py +++ b/setup.py @@ -65,9 +65,7 @@ def run(self): ]) macros = [ - ("USE_COEFFICIENTS", 1), ("NDEBUG", 1), - ("ASSEMBLE_REDUCTION_MATRIX", 1) ] # Robinhood From a661aef03857fdfdb824d97799f77130596e51a9 Mon Sep 17 00:00:00 2001 From: julian Date: Mon, 1 Mar 2021 16:15:25 +0100 Subject: [PATCH 4/5] Fix seutp-minicondav1 issue after github actions update Signed-off-by: julian --- .github/workflows/python-app.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index 4507ca2..f07bd18 100644 --- a/.github/workflows/python-app.yml +++ b/.github/workflows/python-app.yml @@ -51,7 +51,7 @@ jobs: steps: - uses: actions/checkout@v2 - - uses: conda-incubator/setup-miniconda@v1 + - uses: conda-incubator/setup-miniconda@v2 with: auto-update-conda: true python-version: ${{ matrix.python-version }} @@ -76,7 +76,7 @@ jobs: steps: - uses: actions/checkout@v2 - - uses: conda-incubator/setup-miniconda@v1 + - uses: conda-incubator/setup-miniconda@v2 with: auto-update-conda: true python-version: ${{ matrix.python-version }} From 0b267eb309fe5e8530cbd42de7a17462bc829296 Mon Sep 17 00:00:00 2001 From: julian Date: Tue, 2 Mar 2021 09:44:27 +0100 Subject: [PATCH 5/5] [C++] Fix check_overflow from binomial_coeff_table Now that USE_COEFFICIENTS is removed, the conditions in check_overflow needs to support the different check depending on modulus value Signed-off-by: julian --- ripser/ripser.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/ripser/ripser.cpp b/ripser/ripser.cpp index e359a47..d0d16a2 100644 --- a/ripser/ripser.cpp +++ b/ripser/ripser.cpp @@ -83,9 +83,9 @@ static const size_t num_coefficient_bits = 8; static const index_t max_simplex_index = (uintptr_t(1) << (8 * sizeof(index_t) - 1 - num_coefficient_bits)) - 1; -void check_overflow(index_t i) +void check_overflow(const index_t i, const bool is_modulus_2) { - if (i > max_simplex_index) + if ((is_modulus_2 && i < 0) || (!is_modulus_2 && i > max_simplex_index)) throw std::overflow_error( "simplex index " + std::to_string((uint64_t) i) + " in filtration is larger than maximum index " + @@ -99,7 +99,8 @@ class binomial_coeff_table size_t offset; public: - binomial_coeff_table(index_t n, index_t k) : B((n + 1) * (k + 1)) + binomial_coeff_table(const index_t n, const index_t k, + const bool is_modulus_2 = false) : B((n + 1) * (k + 1)) { offset = k + 1; for (index_t i = 0; i <= n; ++i) { @@ -109,7 +110,7 @@ class binomial_coeff_table B[(i - 1) * offset + j - 1] + B[(i - 1) * offset + j]; if (i <= k) B[i * offset + i] = 1; - check_overflow(B[i * offset + std::min(i >> 1, k)]); + check_overflow(B[i * offset + std::min(i >> 1, k)], is_modulus_2); } } @@ -607,7 +608,8 @@ class ripser float _ratio, coefficient_t _modulus, int _do_cocycles) : dist(std::move(_dist)), n(dist.size()), dim_max(_dim_max), threshold(_threshold), ratio(_ratio), modulus(_modulus), - binomial_coeff(n, dim_max + 2), + binomial_coeff(n, dim_max + 2, + std::is_same::value), multiplicative_inverse(multiplicative_inverse_vector(_modulus)), do_cocycles(_do_cocycles) {