From bbde770f2a1a028add4a840ea7291d524081fb83 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Fri, 23 Oct 2015 14:55:18 +0200 Subject: support for finite field coefficients --- ripser.cpp | 94 +++++++++++++++++++++++++++++++++++++++----------------------- 1 file changed, 60 insertions(+), 34 deletions(-) (limited to 'ripser.cpp') diff --git a/ripser.cpp b/ripser.cpp index 75f025d..c8c12ea 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -17,18 +17,19 @@ typedef long index_t; #define USE_BINARY_SEARCH #define USE_EXPONENTIAL_SEARCH -#define INDICATE_PROGRESS +//#define INDICATE_PROGRESS #define PRINT_PERSISTENCE_PAIRS -#define USE_COEFFICIENTS +//#define USE_COEFFICIENTS //#define ASSEMBLE_REDUCTION_COLUMN //#define FILE_FORMAT_DIPHA -#define FILE_FORMAT_UPPER_TRIANGULAR_CSV -//#define FILE_FORMAT_LOWER_TRIANGULAR_CSV +//#define FILE_FORMAT_UPPER_TRIANGULAR_CSV +#define FILE_FORMAT_LOWER_TRIANGULAR_CSV +static const index_t modulus = 2; class binomial_coeff_table { std::vector > B; @@ -169,7 +170,7 @@ struct entry_t { index(_index), value(_value) {} entry_t(index_t _index) : - index(_index), value(0) {} + index(_index), value(1) {} }; entry_t make_entry(index_t _index, index_t _value) { @@ -180,6 +181,10 @@ index_t get_index(entry_t e) { return e.index; } +index_t get_coefficient(entry_t e) { + return e.value; +} + #else typedef index_t entry_t; @@ -188,6 +193,10 @@ index_t get_index(index_t i) { return i; } +index_t get_coefficient(index_t i) { + return 1; +} + entry_t make_entry(index_t _index, index_t _value) { return entry_t(_index); } @@ -446,21 +455,24 @@ public: template inline entry_t pop_pivot(Heap& column) { + #ifndef USE_COEFFICIENTS + static_assert(modulus == 2, "Modulus must be 2 when coefficients are disabled"); + #endif + if( column.empty() ) return entry_t(-1); else { - entry_t max_element = column.top(); - column.pop(); - while( !column.empty() && get_index(column.top()) == get_index(max_element) ) { + index_t pivot_index = get_index(column.top()); + index_t coefficient = 0; + while( !column.empty() && get_index(column.top()) == pivot_index ) { + coefficient = (coefficient + get_coefficient(column.top())) % modulus; + column.pop(); - if( column.empty() ) - return entry_t(-1); - else { - max_element = column.top(); - column.pop(); - } + + if( coefficient == 0 ) + pivot_index = column.empty() ? -1 : get_index(column.top()); } - return max_element; + return make_entry(pivot_index, coefficient); } } @@ -509,7 +521,7 @@ void compute_pairs( std::cout << "persistence intervals in dim " << dim << ":" << std::endl; for (index_t i = 0; i < columns_to_reduce.size(); ++i) { - index_t index = columns_to_reduce[i]; + index_t column_to_reduce = columns_to_reduce[i]; #ifdef ASSEMBLE_REDUCTION_COLUMN std::priority_queue, decltype(comp) > reduction_column(comp); @@ -519,11 +531,14 @@ void compute_pairs( #ifdef INDICATE_PROGRESS std::cout << "\033[K" << "reducing column " << i + 1 << "/" << columns_to_reduce.size() - << " (diameter " << comp_prev.diameter(index) << ")" + << " (diameter " << comp_prev.diameter(column_to_reduce) << ")" << std::flush << "\r"; #endif - index_t pivot, column = index; + index_t column_to_add = column_to_reduce; + + entry_t pivot(column_to_reduce); + std::vector coboundary; @@ -533,47 +548,56 @@ void compute_pairs( reduction_column.push( column ); #endif - simplex_coboundary_enumerator coboundaries(column, dim, n, binomial_coeff); + simplex_coboundary_enumerator coboundaries(column_to_add, dim, n, binomial_coeff); while (coboundaries.has_next()) { entry_t coface = coboundaries.next(); index_t coface_index = get_index(coface); if (comp.diameter(coface_index) <= threshold) - working_coboundary.push(coface); + working_coboundary.push(make_entry(coface_index, get_coefficient(pivot))); } - pivot = get_index(get_pivot(working_coboundary)); + pivot = get_pivot(working_coboundary); - if (pivot != -1) { - auto pair = pivot_column_index.find(pivot); + if (get_index(pivot) != -1) { + auto pair = pivot_column_index.find(get_index(pivot)); if (pair == pivot_column_index.end()) { - pivot_column_index.insert(std::make_pair(pivot, index)); + pivot_column_index.insert(std::make_pair(get_index(pivot), column_to_reduce)); #ifdef PRINT_PERSISTENCE_PAIRS - value_t birth = comp_prev.diameter(index), death = comp.diameter(pivot); - if (birth != death) - std::cout << "\033[K" << " [" << birth << "," << death << ")" << std::endl << std::flush; + value_t birth = comp_prev.diameter(column_to_reduce), death = comp.diameter(get_index(pivot)); + if (birth != death) { + #ifdef INDICATE_PROGRESS + std::cout << "\033[K"; + #endif + std::cout << " [" << birth << "," << death << ")" << std::endl << std::flush; + } #endif break; } - column = pair->second; + column_to_add = pair->second; } - } while ( pivot != -1 ); + } while ( get_index(pivot) != -1 ); #ifdef PRINT_PERSISTENCE_PAIRS - if ( pivot == -1 ) { - value_t birth = comp_prev.diameter(index); - std::cout << "\033[K" << " [" << birth << ", )" << std::endl << std::flush; + if ( get_index(pivot) == -1 ) { + value_t birth = comp_prev.diameter(column_to_reduce); + #ifdef INDICATE_PROGRESS + std::cout << "\033[K"; + #endif + std::cout << " [" << birth << ", )" << std::endl << std::flush; } #endif } + #ifdef INDICATE_PROGRESS std::cout << "\033[K"; + #endif } void print_usage_and_exit(int exit_code) @@ -692,12 +716,12 @@ int main( int argc, char** argv ) { #ifdef FILE_FORMAT_LOWER_TRIANGULAR_CSV - std::vector distances; + std::vector& distances = diameters[1]; std::string value_string; while(std::getline(input_stream, value_string, ',')) distances.push_back(std::stod(value_string)); - compressed_lower_distance_matrix dist(std::move(distances)); + compressed_lower_distance_matrix_adapter dist(distances); index_t n = dist.size(); @@ -705,6 +729,8 @@ int main( int argc, char** argv ) { #endif + + auto result = std::minmax_element(dist.distances.begin(), dist.distances.end()); std::cout << "value range: [" << *result.first << "," << *result.second << "]" << std::endl; -- cgit v1.2.3