diff options
-rw-r--r-- | Makefile | 7 | ||||
-rw-r--r-- | ripser.cpp | 124 |
2 files changed, 68 insertions, 63 deletions
@@ -1,15 +1,18 @@ build: ripser -all: ripser ripser-debug +all: ripser ripser-coeff ripser-debug ripser: ripser.cpp c++ -std=c++11 ripser.cpp -o ripser -Ofast -D NDEBUG +ripser-coeff: ripser.cpp + c++ -std=c++11 ripser.cpp -o ripser-coeff -Ofast -D NDEBUG -D USE_COEFFICIENTS + ripser-debug: ripser.cpp c++ -std=c++11 ripser.cpp -o ripser-debug -g clean: - rm -f ripser ripser-debug + rm -f ripser ripser-coeff ripser-debug @@ -38,7 +38,7 @@ //#define USE_COEFFICIENTS -#define INDICATE_PROGRESS +//#define INDICATE_PROGRESS #define PRINT_PERSISTENCE_PAIRS //#define USE_GOOGLE_HASHMAP @@ -73,8 +73,10 @@ static const std::chrono::milliseconds time_step(40); static const std::string clear_line("\r\033[K"); +static const size_t num_coefficient_bits = 8; + static const index_t max_simplex_index = - (1l << (8 * (sizeof(index_t) - sizeof(coefficient_t)) - 1)) - 1; + (1l << (8 * sizeof(index_t) - 1 - num_coefficient_bits)) - 1; void check_overflow(index_t i) { if @@ -84,7 +86,7 @@ void check_overflow(index_t i) { (i < 0) #endif throw std::overflow_error("simplex index " + std::to_string((uint64_t)i) + - " in filtration is larger than maximum index" + + " in filtration is larger than maximum index " + std::to_string(max_simplex_index)); } @@ -129,8 +131,8 @@ std::vector<coefficient_t> multiplicative_inverse_vector(const coefficient_t m) #ifdef USE_COEFFICIENTS struct __attribute__((packed)) entry_t { - index_t index : 8 * (sizeof(index_t) - sizeof(coefficient_t)); - coefficient_t coefficient; + index_t index : 8 * sizeof(index_t) - num_coefficient_bits; + coefficient_t coefficient : num_coefficient_bits; entry_t(index_t _index, coefficient_t _coefficient) : index(_index), coefficient(_coefficient) {} entry_t(index_t _index) : index(_index), coefficient(0) {} @@ -144,9 +146,6 @@ 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; } -bool operator==(const entry_t& e1, const entry_t& e2) { - return get_index(e1) == get_index(e2) && get_coefficient(e1) == get_coefficient(e2); -} std::ostream& operator<<(std::ostream& stream, const entry_t& e) { stream << get_index(e) << ":" << get_coefficient(e); @@ -165,18 +164,6 @@ void set_coefficient(entry_t& e, const coefficient_t c) {} const entry_t& get_entry(const entry_t& e) { return e; } -struct entry_hash { - std::size_t operator()(const entry_t& e) const { return std::hash<index_t>()(get_index(e)); } -}; - -struct equal_index { - bool operator()(const entry_t& e, const entry_t& f) const { - return get_index(e) == get_index(f); - } -}; - -typedef hash_map<entry_t, index_t, entry_hash, equal_index> entry_hash_map; - typedef std::pair<value_t, index_t> diameter_index_t; value_t get_diameter(const diameter_index_t& i) { return i.first; } index_t get_index(const diameter_index_t& i) { return i.second; } @@ -344,41 +331,6 @@ public: } }; -template <typename Column> diameter_entry_t pop_pivot(Column& column, coefficient_t modulus) { - diameter_entry_t pivot(-1); -#ifdef USE_COEFFICIENTS - while (!column.empty()) { - if (get_coefficient(pivot) == 0) - pivot = column.top(); - else if (get_index(column.top()) != get_index(pivot)) - break; - else - set_coefficient(pivot, - (get_coefficient(pivot) + get_coefficient(column.top())) % modulus); - column.pop(); - } - if (get_coefficient(pivot) == 0) pivot = -1; -#else - while (!column.empty()) { - pivot = column.top(); - column.pop(); - if (!column.empty()) { - if (get_index(column.top()) != get_index(pivot)) - return pivot; - column.pop(); - } - } - pivot = -1; -#endif - return pivot; -} - -template <typename Column> diameter_entry_t get_pivot(Column& column, const coefficient_t modulus) { - diameter_entry_t result = pop_pivot(column, modulus); - if (get_index(result) != -1) column.push(result); - return result; -} - template <typename T> T begin(std::pair<T, T>& p) { return p.first; } template <typename T> T end(std::pair<T, T>& p) { return p.second; } @@ -431,6 +383,19 @@ template <typename DistanceMatrix> class ripser { const binomial_coeff_table binomial_coeff; const std::vector<coefficient_t> multiplicative_inverse; mutable std::vector<diameter_entry_t> coface_entries; + + struct entry_hash { + std::size_t operator()(const entry_t& e) const { return std::hash<index_t>()(::get_index(e)); } + }; + + struct equal_index { + bool operator()(const entry_t& e, const entry_t& f) const { + return ::get_index(e) == ::get_index(f); + } + }; + + typedef hash_map<entry_t, index_t, entry_hash, equal_index> entry_hash_map; + public: ripser(DistanceMatrix&& _dist, index_t _dim_max, value_t _threshold, float _ratio, @@ -488,11 +453,13 @@ public: } #endif auto coface = cofaces.next(); + if (get_diameter(coface) <= threshold) { - next_simplices.push_back({get_diameter(coface), get_index(coface)}); + next_simplices.push_back({get_diameter(coface), get_index(coface)}); - if (pivot_column_index.find(get_entry(coface)) == pivot_column_index.end()) - columns_to_reduce.push_back({get_diameter(coface), get_index(coface)}); + if (pivot_column_index.find(get_entry(coface)) == pivot_column_index.end()) + columns_to_reduce.push_back({get_diameter(coface), get_index(coface)}); + } } } @@ -543,6 +510,40 @@ public: #endif } + template <typename Column> diameter_entry_t pop_pivot(Column& column) { + diameter_entry_t pivot(-1); +#ifdef USE_COEFFICIENTS + while (!column.empty()) { + if (get_coefficient(pivot) == 0) + pivot = column.top(); + else if (get_index(column.top()) != get_index(pivot)) + break; + else + set_coefficient(pivot, + (get_coefficient(pivot) + get_coefficient(column.top())) % modulus); + column.pop(); + } + if (get_coefficient(pivot) == 0) pivot = -1; +#else + while (!column.empty()) { + pivot = column.top(); + column.pop(); + if (!column.empty()) { + if (get_index(column.top()) != get_index(pivot)) return pivot; + column.pop(); + } + } + pivot = -1; +#endif + return pivot; + } + + template <typename Column> diameter_entry_t get_pivot(Column& column) { + diameter_entry_t result = pop_pivot(column); + if (get_index(result) != -1) column.push(result); + return result; + } + template <typename Column> diameter_entry_t init_coboundary_and_get_pivot(const diameter_entry_t simplex, Column& working_coboundary, const index_t& dim, @@ -562,7 +563,7 @@ public: } } for (auto coface : coface_entries) working_coboundary.push(coface); - return get_pivot(working_coboundary, modulus); + return get_pivot(working_coboundary); } template <typename Column> @@ -640,7 +641,7 @@ public: add_coboundary(reduction_matrix, columns_to_reduce, index_column_to_add, factor, working_reduction_column, working_coboundary, dim); - pivot = get_pivot(working_coboundary, modulus); + pivot = get_pivot(working_coboundary); } else { #ifdef PRINT_PERSISTENCE_PAIRS value_t death = get_diameter(pivot); @@ -654,7 +655,7 @@ public: pivot_column_index.insert({get_entry(pivot), index_column_to_reduce}); while (true) { - diameter_entry_t e = pop_pivot(working_reduction_column, modulus); + diameter_entry_t e = pop_pivot(working_reduction_column); if (get_index(e) == -1) break; assert(get_coefficient(e) > 0); reduction_matrix.push_back(e); @@ -1123,5 +1124,6 @@ int main(int argc, char** argv) { dim_max, threshold, ratio, modulus) .compute_barcodes(); } + exit(0); } } |