From 446187d05480fddcb6de2919b8ac10c6cc19212c Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 15 Nov 2015 17:10:03 -0500 Subject: build macros in project settings --- ripser.xcodeproj/project.pbxproj | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'ripser.xcodeproj/project.pbxproj') diff --git a/ripser.xcodeproj/project.pbxproj b/ripser.xcodeproj/project.pbxproj index 196bac0..9a3f1e5 100644 --- a/ripser.xcodeproj/project.pbxproj +++ b/ripser.xcodeproj/project.pbxproj @@ -205,6 +205,10 @@ 551018501BD63CB300990BFF /* Debug */ = { isa = XCBuildConfiguration; buildSettings = { + GCC_PREPROCESSOR_DEFINITIONS = ( + "$(inherited)", + FILE_FORMAT_LOWER_TRIANGULAR_CSV, + ); PRODUCT_NAME = "$(TARGET_NAME)"; }; name = Debug; -- cgit v1.2.3 From 6698973093f75721d2e188da279fd799e753c4e9 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Sun, 15 Nov 2015 18:00:15 -0500 Subject: coface diameter speedup --- ripser.cpp | 39 +++++++++++++++++++++++++++++++-------- ripser.xcodeproj/project.pbxproj | 1 + 2 files changed, 32 insertions(+), 8 deletions(-) (limited to 'ripser.xcodeproj/project.pbxproj') diff --git a/ripser.cpp b/ripser.cpp index 0ed66e8..902caed 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -321,7 +321,7 @@ public: return k != -1 && n != -1; } - inline entry_t next() + inline std::pair next() { while ( binomial_coeff( n , k ) <= idx ) { idx -= binomial_coeff( n , k ); @@ -331,10 +331,12 @@ public: --n; } - return make_entry( - modified_idx + binomial_coeff( n-- , k + 1 ), + auto result = std::make_pair(make_entry( + modified_idx + binomial_coeff( n , k + 1 ), k & 1 ? -1 : 1 - ); + ), n); + --n; + return result; } }; @@ -360,7 +362,7 @@ std::vector get_diameters ( simplex_coboundary_enumerator coboundaries(simplex, dim - 1, n, binomial_coeff); while (coboundaries.has_next()) { - index_t coface = get_index(coboundaries.next()); + index_t coface = get_index(coboundaries.next().first); diameters[coface] = std::max( diameters[coface], previous_diameters[simplex]); } } @@ -750,10 +752,11 @@ inline void push_entry(Heap& column, index_t i, coefficient_t c, value_t diamete } -template +template void compute_pairs( std::vector& columns_to_reduce, std::unordered_map& pivot_column_index, + const DistanceMatrix& dist, const ComparatorCofaces& comp, const Comparator& comp_prev, index_t dim, index_t n, value_t threshold, @@ -863,12 +866,30 @@ void compute_pairs( #endif #endif + //TODO: cache instead of compute? + value_t simplex_diameter = comp_prev.diameter(get_index(simplex)); + + + std::vector vertices = vertices_of_simplex(get_index(simplex), dim, n, binomial_coeff); + simplex_coboundary_enumerator cofaces(get_index(simplex), dim, n, binomial_coeff); while (cofaces.has_next()) { - entry_t coface = cofaces.next(); + auto coface_descriptor = cofaces.next(); + entry_t coface = coface_descriptor.first; + index_t covertex = coface_descriptor.second; +// index_t coface_index = get_index(coface); - value_t coface_diameter = comp.diameter(coface_index); +// + value_t coface_diameter = simplex_diameter; + for (index_t v : vertices) { + coface_diameter = std::max(coface_diameter, dist(v, covertex)); + } + +// value_t coface_diameter = comp.diameter(coface_index); + +// assert(coface_diameter_variant == coface_diameter); + if (coface_diameter <= threshold) { coefficient_t coface_coefficient = get_coefficient(coface) + modulus; coface_coefficient %= modulus; @@ -1181,6 +1202,7 @@ int main( int argc, char** argv ) { compute_pairs( columns_to_reduce, pivot_column_index, + dist, comp, comp_prev, dim, n, threshold, @@ -1231,6 +1253,7 @@ int main( int argc, char** argv ) { compute_pairs( columns_to_reduce, pivot_column_index, + dist, comp, comp_prev, dim, n, threshold, diff --git a/ripser.xcodeproj/project.pbxproj b/ripser.xcodeproj/project.pbxproj index 9a3f1e5..430ef6b 100644 --- a/ripser.xcodeproj/project.pbxproj +++ b/ripser.xcodeproj/project.pbxproj @@ -208,6 +208,7 @@ GCC_PREPROCESSOR_DEFINITIONS = ( "$(inherited)", FILE_FORMAT_LOWER_TRIANGULAR_CSV, + STORE_DIAMETERS, ); PRODUCT_NAME = "$(TARGET_NAME)"; }; -- cgit v1.2.3 From 6812ded36b5ada39c7f16a67f3235367a3c67f4a Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Mon, 16 Nov 2015 00:48:46 -0500 Subject: store diameters in reduction matrix --- ripser.cpp | 140 +++++++++++++++++++++------------------ ripser.xcodeproj/project.pbxproj | 2 + 2 files changed, 76 insertions(+), 66 deletions(-) (limited to 'ripser.xcodeproj/project.pbxproj') diff --git a/ripser.cpp b/ripser.cpp index 902caed..7fbbe66 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -193,43 +193,62 @@ inline void set_coefficient(index_t& e, const coefficient_t c) { e = c; } inline const entry_t& get_entry(const entry_t& e) { return e; } +template struct greater_index { - bool operator() (const entry_t& a, const entry_t& b) { + bool operator() (const Entry& a, const Entry& b) { return get_index(a) > get_index(b); } }; #ifdef STORE_DIAMETERS -class entry_diameter_t: public std::pair { +typedef std::pair diameter_index_t; +inline value_t get_diameter(diameter_index_t i) { + return i.first; +} +inline index_t get_index(diameter_index_t i) { + return i.second; +} +#else +typedef index_t diameter_index_t; +#endif + +#ifdef STORE_DIAMETERS +class diameter_entry_t: public std::pair { public: - entry_diameter_t(std::pair p) : std::pair(p) {} + diameter_entry_t(std::pair p) : std::pair(p) {} #ifdef USE_COEFFICIENTS - entry_diameter_t(entry_t e) : std::pair(e, 0) {} + diameter_entry_t(entry_t e) : std::pair(0, e) {} #endif - entry_diameter_t(index_t i) : std::pair(i, 0) {} + diameter_entry_t(index_t i) : std::pair(0, i) {} }; -inline const entry_t& get_entry(const entry_diameter_t& p) { return p.first; } -inline entry_t& get_entry(entry_diameter_t& p) { return p.first; } -inline const index_t get_index(const entry_diameter_t& p) { return get_index(get_entry(p)); } -inline const coefficient_t get_coefficient(const entry_diameter_t& p) { return get_coefficient(get_entry(p)); } -inline const value_t& get_diameter(const entry_diameter_t& p) { return p.second; } -inline void set_coefficient(entry_diameter_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); } -inline entry_diameter_t make_entry_diameter(index_t _index, coefficient_t _coefficient, value_t _diameter) { - return std::make_pair(make_entry(_index, _coefficient), _diameter); +inline const entry_t& get_entry(const diameter_entry_t& p) { return p.second; } +inline entry_t& get_entry(diameter_entry_t& p) { return p.second; } +inline const index_t get_index(const diameter_entry_t& p) { return get_index(get_entry(p)); } +inline const coefficient_t get_coefficient(const diameter_entry_t& p) { return get_coefficient(get_entry(p)); } +inline const value_t& get_diameter(const diameter_entry_t& p) { return p.first; } +inline void set_coefficient(diameter_entry_t& p, const coefficient_t c) { set_coefficient(get_entry(p), c); } +inline diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index, coefficient_t _coefficient) { + return std::make_pair(_diameter, make_entry(_index, _coefficient)); +} +inline diameter_entry_t make_diameter_entry(diameter_index_t _diameter_index, coefficient_t _coefficient) { + return std::make_pair(get_diameter(_diameter_index), make_entry(get_index(_diameter_index), _coefficient)); } struct greater_diameter_or_index { - inline bool operator() (const entry_diameter_t& a, const entry_diameter_t& b) { + inline bool operator() (const diameter_entry_t& a, const diameter_entry_t& b) { return (get_diameter(a) > get_diameter(b)) || ((get_diameter(a) == get_diameter(b)) && (get_index(a) > get_index(b))); } }; #else -typedef entry_t entry_diameter_t; -inline entry_diameter_t make_entry_diameter(index_t _index, coefficient_t _coefficient, value_t _diameter) { +typedef entry_t diameter_entry_t; +inline diameter_entry_t make_diameter_entry(value_t _diameter, index_t _index, coefficient_t _coefficient) { + return make_entry(_index, _coefficient); +} +inline diameter_entry_t make_diameter_entry(index_t _index, coefficient_t _coefficient) { return make_entry(_index, _coefficient); } #endif @@ -333,8 +352,9 @@ public: } auto result = std::make_pair(make_entry( modified_idx + binomial_coeff( n , k + 1 ), - k & 1 ? -1 : 1 - ), n); + k & 1 ? -1 : 1), + n); + --n; return result; } @@ -552,7 +572,7 @@ public: #ifdef USE_COEFFICIENTS template -inline entry_diameter_t pop_pivot(Heap& column, coefficient_t modulus) +inline diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) { if( column.empty() ) @@ -583,7 +603,7 @@ inline entry_diameter_t pop_pivot(Heap& column, coefficient_t modulus) #else template -inline entry_diameter_t pop_pivot(Heap& column, coefficient_t modulus) +inline diameter_entry_t pop_pivot(Heap& column, coefficient_t modulus) { if( column.empty() ) @@ -607,9 +627,9 @@ inline entry_diameter_t pop_pivot(Heap& column, coefficient_t modulus) #endif template -inline entry_diameter_t get_pivot(Heap& column, coefficient_t modulus) +inline diameter_entry_t get_pivot(Heap& column, coefficient_t modulus) { - entry_diameter_t result = pop_pivot(column, modulus); + diameter_entry_t result = pop_pivot(column, modulus); if( get_index(result) != -1 ) { column.push(result); } @@ -617,17 +637,6 @@ inline entry_diameter_t get_pivot(Heap& column, coefficient_t modulus) } -#ifdef STORE_DIAMETERS -typedef std::pair diameter_index_t; -inline value_t get_diameter(diameter_index_t i) { - return i.first; -} -inline index_t get_index(diameter_index_t i) { - return i.second; -} -#else -typedef index_t diameter_index_t; -#endif template void assemble_columns_to_reduce ( @@ -745,7 +754,7 @@ template inline void push_entry(Heap& column, index_t i, coefficient_t c, value_t diameter) { entry_t e = make_entry(i, c); #ifdef STORE_DIAMETERS - column.push(std::make_pair(e, diameter)); + column.push(std::make_pair(diameter, e)); #else column.push(e); #endif @@ -769,27 +778,27 @@ void compute_pairs( #endif #ifdef ASSEMBLE_REDUCTION_MATRIX - compressed_sparse_matrix reduction_matrix; + compressed_sparse_matrix reduction_matrix; #else #ifdef USE_COEFFICIENTS - std::vector reduction_coefficients; + std::vector reduction_coefficients; #endif #endif // size_t boundary_additions = 0; for (index_t i = 0; i < columns_to_reduce.size(); ++i) { - index_t column_to_reduce = get_index(columns_to_reduce[i]); + auto column_to_reduce = columns_to_reduce[i]; #ifdef ASSEMBLE_REDUCTION_MATRIX - std::priority_queue, greater_index> reduction_column; + std::priority_queue, greater_index> reduction_column; #endif #ifdef STORE_DIAMETERS - std::priority_queue, greater_diameter_or_index > working_coboundary; + std::priority_queue, greater_diameter_or_index > working_coboundary; #else std::priority_queue, decltype(comp) > working_coboundary(comp); #endif @@ -797,7 +806,7 @@ void compute_pairs( #ifdef STORE_DIAMETERS value_t diameter = get_diameter(columns_to_reduce[i]); #else - value_t diameter = comp_prev.diameter(column_to_reduce); + value_t diameter = comp_prev.diameter(get_index(column_to_reduce)); #endif @@ -809,13 +818,13 @@ void compute_pairs( index_t j = i; - index_t column_to_add = column_to_reduce; + auto column_to_add = column_to_reduce; // start with a dummy pivot entry with coefficient -1 in order to initialize working_coboundary // with the coboundary of the simplex with index column_to_reduce - entry_diameter_t pivot = make_entry_diameter(-1, modulus - 1, 0); + diameter_entry_t pivot = make_diameter_entry(0, -1, modulus - 1); // std::cout << "reducing " << column_to_reduce << ":" << std::endl; @@ -826,10 +835,10 @@ void compute_pairs( // initialize reduction_matrix as identity matrix #ifdef ASSEMBLE_REDUCTION_MATRIX - reduction_matrix.push_back(make_entry(column_to_reduce, 1)); + reduction_matrix.push_back(make_diameter_entry(column_to_reduce, 1)); #else #ifdef USE_COEFFICIENTS - reduction_coefficients.push_back(make_entry(column_to_reduce, 1)); + reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, 1)); #endif #endif @@ -840,7 +849,7 @@ void compute_pairs( const coefficient_t factor = modulus - get_coefficient(pivot); -// std::priority_queue, decltype(comp) > eliminating_coboundary(comp); +// std::priority_queue, decltype(comp) > eliminating_coboundary(comp); // std::cout << "w:" << get_column_vector(working_coboundary, modulus) << std::endl; @@ -853,22 +862,28 @@ void compute_pairs( #ifdef ASSEMBLE_REDUCTION_MATRIX - entry_t simplex = *it; + auto simplex = *it; coefficient_t simplex_coefficient = get_coefficient(simplex); simplex_coefficient *= factor; simplex_coefficient %= modulus; - reduction_column.push( make_entry( get_index(simplex), simplex_coefficient ) ); #else #ifdef USE_COEFFICIENTS - const entry_t& simplex = reduction_coefficients[j]; + const auto& simplex = reduction_coefficients[j]; #else - const entry_t simplex = column_to_add; + const auto simplex = column_to_add; #endif #endif - //TODO: cache instead of compute? + #ifdef STORE_DIAMETERS + value_t simplex_diameter = get_diameter(simplex); + #else value_t simplex_diameter = comp_prev.diameter(get_index(simplex)); + #endif + assert(simplex_diameter == comp_prev.diameter(get_index(simplex))); + #ifdef ASSEMBLE_REDUCTION_MATRIX + reduction_column.push( make_diameter_entry( simplex_diameter, get_index(simplex), simplex_coefficient ) ); + #endif std::vector vertices = vertices_of_simplex(get_index(simplex), dim, n, binomial_coeff); @@ -878,27 +893,20 @@ void compute_pairs( entry_t coface = coface_descriptor.first; index_t covertex = coface_descriptor.second; -// + index_t coface_index = get_index(coface); -// + value_t coface_diameter = simplex_diameter; for (index_t v : vertices) { coface_diameter = std::max(coface_diameter, dist(v, covertex)); } -// value_t coface_diameter = comp.diameter(coface_index); - -// assert(coface_diameter_variant == coface_diameter); + assert(comp.diameter(coface_index) == coface_diameter); if (coface_diameter <= threshold) { - coefficient_t coface_coefficient = get_coefficient(coface) + modulus; - coface_coefficient %= modulus; - - coface_coefficient *= get_coefficient(simplex); - coface_coefficient %= modulus; - - coface_coefficient *= factor; + coefficient_t coface_coefficient = get_coefficient(coface) * get_coefficient(simplex) * factor; coface_coefficient %= modulus; + if (coface_coefficient < 0) coface_coefficient += modulus; assert(coface_coefficient >= 0); @@ -941,18 +949,18 @@ void compute_pairs( #ifdef ASSEMBLE_REDUCTION_MATRIX reduction_matrix.pop_back(); while (true) { - entry_t e = get_entry(pop_pivot(reduction_column, modulus)); + diameter_entry_t e = pop_pivot(reduction_column, modulus); index_t index = get_index(e); if (index == -1) break; const coefficient_t coefficient = inverse * get_coefficient(e) % modulus; assert(coefficient > 0); - reduction_matrix.push_back(make_entry(index, coefficient)); + reduction_matrix.push_back(make_diameter_entry(get_diameter(e), index, coefficient)); } #else #ifdef USE_COEFFICIENTS reduction_coefficients.pop_back(); - reduction_coefficients.push_back(make_entry(column_to_reduce, inverse)); + reduction_coefficients.push_back(make_diameter_entry(column_to_reduce, inverse)); #endif #endif @@ -977,7 +985,7 @@ void compute_pairs( } j = pair->second; - column_to_add = get_index(columns_to_reduce[j]); + column_to_add = columns_to_reduce[j]; } } while ( get_index(pivot) != -1 ); diff --git a/ripser.xcodeproj/project.pbxproj b/ripser.xcodeproj/project.pbxproj index 430ef6b..290272d 100644 --- a/ripser.xcodeproj/project.pbxproj +++ b/ripser.xcodeproj/project.pbxproj @@ -208,7 +208,9 @@ GCC_PREPROCESSOR_DEFINITIONS = ( "$(inherited)", FILE_FORMAT_LOWER_TRIANGULAR_CSV, + USE_COEFFICIENTS, STORE_DIAMETERS, + PRINT_PERSISTENCE_PAIRS, ); PRODUCT_NAME = "$(TARGET_NAME)"; }; -- cgit v1.2.3