From c86ac6dd00ee59c096646cc2e192d6daf529e631 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 22 Oct 2015 10:50:52 +0200 Subject: binary search for vertex enumeration --- ripser.cpp | 109 +++++++++++++++++++++++++++++++++++++++++-------------------- 1 file changed, 73 insertions(+), 36 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 95ea159..137d07b 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -14,7 +14,13 @@ typedef long index_t; #define PRECOMPUTE_DIAMETERS -#define PRECOMPUTE_DIAMETERS_IN_TOP_DIMENSION +//#define PRECOMPUTE_DIAMETERS_IN_TOP_DIMENSION + +#define USE_BINARY_SEARCH + +#define INDICATE_PROGRESS + +#define PRINT_PERSISTENCE_PAIRS //#define ASSEMBLE_REDUCTION_COLUMN @@ -46,7 +52,7 @@ public: } } - index_t operator()(index_t n, index_t k) const { + inline index_t operator()(index_t n, index_t k) const { // std::cout << "B(" << n << "," << k << ")\n"; assert(n <= n_max); assert(k <= k_max); @@ -57,14 +63,36 @@ public: template -OutputIterator get_simplex_vertices( index_t idx, index_t dim, index_t n, const binomial_coeff_table& binomial_coeff, OutputIterator out ) +OutputIterator get_simplex_vertices( index_t idx, const index_t dim, index_t n, const binomial_coeff_table& binomial_coeff, OutputIterator out ) { --n; for( index_t k = dim + 1; k > 0; --k ) { - while( binomial_coeff( n , k ) > idx ) { - --n; + + #ifdef USE_BINARY_SEARCH + if ( binomial_coeff( n , k ) > idx ) { + index_t count; + + for (count = 1; (binomial_coeff( n - count , k ) > idx); count = std::min(count << 1, n)); + + while (count > 0) { + index_t i = n; + index_t step = count >> 1; + i -= step; + if (binomial_coeff( i , k ) > idx) { + n = --i; + count -= step + 1; + } else count = step; + } } + #else + while( binomial_coeff( n , k ) > idx ) + --n; + #endif + + assert( binomial_coeff( n , k ) <= idx ); + assert( binomial_coeff( n + 1, k ) > idx ); + *out++ = n; idx -= binomial_coeff( n , k ); @@ -74,7 +102,7 @@ OutputIterator get_simplex_vertices( index_t idx, index_t dim, index_t n, const return out; } -std::vector vertices_of_simplex(index_t simplex_index, index_t dim, index_t n, const binomial_coeff_table& binomial_coeff) { +std::vector vertices_of_simplex(const index_t simplex_index, const index_t dim, const index_t n, const binomial_coeff_table& binomial_coeff) { std::vector vertices; get_simplex_vertices( simplex_index, dim, n, binomial_coeff, std::back_inserter(vertices) ); return vertices; @@ -110,7 +138,7 @@ public: // reverse(_reverse) {}; - dist_t diameter(const index_t index) const { + inline dist_t diameter(const index_t index) const { dist_t diam = 0; get_simplex_vertices(index, dim, dist.size(), binomial_coeff, vertices.begin() ); @@ -121,7 +149,7 @@ public: return diam; } - bool operator()(const index_t a, const index_t b) const + inline bool operator()(const index_t a, const index_t b) const { assert(a < binomial_coeff(dist.size(), dim + 1)); assert(b < binomial_coeff(dist.size(), dim + 1)); @@ -136,13 +164,10 @@ public: for (index_t j = i + 1; j <= dim; ++j) { a_diam = std::max(a_diam, dist(vertices[i], vertices[j])); if (a_diam > b_diam) -// (((a_diam < b_diam) && !reverse) || -// ((a_diam > b_diam) && reverse)) return true; } return (a_diam == b_diam) && (a > b); -// return (a_diam == b_diam) && (((a < b) && !reverse) || ((a > b) && reverse)); } }; @@ -173,7 +198,7 @@ public: binomial_coeff(_binomial_coeff) {}; - dist_t diameter(const index_t index) const { + inline dist_t diameter(const index_t index) const { dist_t diam = 0; get_simplex_vertices(index, 2, dist.size(), binomial_coeff, vertices.begin() ); @@ -187,7 +212,7 @@ public: return diam; } - bool operator()(const index_t a, const index_t b) const + inline bool operator()(const index_t a, const index_t b) const { assert(a < binomial_coeff(dist.size(), 3)); assert(b < binomial_coeff(dist.size(), 3)); @@ -212,7 +237,7 @@ public: template -void get_simplex_coboundary( index_t idx, index_t dim, index_t n, const binomial_coeff_table& binomial_coeff, OutputIterator coboundary ) +void inline get_simplex_coboundary( index_t idx, index_t dim, index_t n, const binomial_coeff_table& binomial_coeff, OutputIterator coboundary ) { --n; @@ -244,7 +269,7 @@ private: const binomial_coeff_table& binomial_coeff; public: - simplex_coboundary_enumerator ( + inline simplex_coboundary_enumerator ( index_t _idx, index_t _dim, index_t _n, @@ -254,7 +279,7 @@ public: } - bool has_next() + inline bool has_next() { while ( (k != -1 && n != -1) && (binomial_coeff( n , k ) <= idx) ) { idx -= binomial_coeff( n , k ); @@ -268,7 +293,7 @@ public: return k != -1 && n != -1; } - index_t next() + inline index_t next() { while ( binomial_coeff( n , k ) <= idx ) { idx -= binomial_coeff( n , k ); @@ -299,8 +324,9 @@ std::vector get_diameters ( for (index_t simplex = 0; simplex < previous_diameters.size(); ++simplex) { coboundary.clear(); + #ifdef INDICATE_PROGRESS std::cout << "\033[Kpropagating diameter of simplex " << simplex + 1 << "/" << previous_diameters.size() << std::flush << "\r"; - + #endif simplex_coboundary_enumerator coboundaries(simplex, dim - 1, n, binomial_coeff); while (coboundaries.has_next()) { @@ -309,7 +335,9 @@ std::vector get_diameters ( } } + #ifdef INDICATE_PROGRESS std::cout << "\033[K"; + #endif // rips_filtration_comparator comp(dist, dim, binomial_coeff); // for (index_t simplex = 0; simplex < diameters.size(); ++simplex) { @@ -344,12 +372,12 @@ public: binomial_coeff(_binomial_coeff) {}; - value_t diameter(const index_t a) const { + inline value_t diameter(const index_t a) const { assert(a < diameters.size()); return diameters[a]; } - bool operator()(const index_t a, const index_t b) const + inline bool operator()(const index_t a, const index_t b) const { assert(a < diameters.size()); assert(b < diameters.size()); @@ -366,11 +394,11 @@ class distance_matrix { public: typedef value_t value_type; std::vector > distances; - value_t operator()(const index_t a, const index_t b) const { + inline value_t operator()(const index_t a, const index_t b) const { return distances[a][b]; } - size_t size() const { + inline size_t size() const { return distances.size(); } @@ -399,7 +427,7 @@ public: } } - value_t operator()(const index_t a, const index_t b) const { + inline value_t operator()(const index_t a, const index_t b) const { if (a < b) return rows[a][b]; else if (a > b) @@ -408,7 +436,7 @@ public: return 0; } - size_t size() const { + inline size_t size() const { return n; } @@ -421,22 +449,22 @@ public: std::vector entries; - size_t size() const { + inline size_t size() const { return bounds.size(); } - typename std::vector::const_iterator cbegin(size_t index) { + inline typename std::vector::const_iterator cbegin(size_t index) const { assert(index < size()); return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1]; } - typename std::vector::const_iterator cend(size_t index) { + inline typename std::vector::const_iterator cend(size_t index) const { assert(index < size()); return entries.cbegin() + bounds[index]; } template - void append(Iterator begin, Iterator end) { + inline void append(Iterator begin, Iterator end) { for (Iterator it = begin; it != end; ++it) { entries.push_back(*it); } @@ -445,14 +473,14 @@ public: template - void append(const Collection collection) { + inline void append(const Collection collection) { append(collection.cbegin(), collection.cend()); } }; template -index_t pop_pivot(Heap& column) +inline index_t pop_pivot(Heap& column) { if( column.empty() ) return -1; @@ -473,7 +501,7 @@ index_t pop_pivot(Heap& column) } template -index_t get_pivot(Heap& column) +inline index_t get_pivot(Heap& column) { index_t max_element = pop_pivot(column); if( max_element != -1 ) @@ -482,7 +510,7 @@ index_t get_pivot(Heap& column) } template -std::vector move_to_column_vector(Heap& column) +inline std::vector move_to_column_vector(Heap& column) { std::vector temp_col; index_t pivot = pop_pivot( column ); @@ -529,9 +557,11 @@ void compute_pairs( std::priority_queue, decltype(comp) > working_coboundary(comp); + #ifdef INDICATE_PROGRESS std::cout << "\033[K" << "reducing column " << i + 1 << "/" << columns_to_reduce.size() << " (diameter " << comp_prev.diameter(index) << ")" << std::flush << "\r"; + #endif index_t pivot, column = index; @@ -567,9 +597,11 @@ void compute_pairs( if (pair == pivot_column_index.end()) { pivot_column_index.insert(std::make_pair(pivot, index)); + #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; + #endif break; } @@ -579,11 +611,12 @@ void compute_pairs( } while ( 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; } + #endif } @@ -623,7 +656,7 @@ int main( int argc, char** argv ) { print_help_and_exit(); std::string parameter = std::string( argv[ idx ] ); size_t pos_last_digit; - threshold = std::stod( parameter, &pos_last_digit ); + threshold = std::stof( parameter, &pos_last_digit ); if( pos_last_digit != parameter.size() ) print_help_and_exit(); } else print_help_and_exit(); @@ -733,14 +766,18 @@ int main( int argc, char** argv ) { std::vector edge_vertices(2); for (index_t edge = 0; edge < diameters.size(); ++edge) { + #ifdef INDICATE_PROGRESS std::cout << "\033[Kstoring diameter of simplex " << edge + 1 << "/" << diameters.size() << std::flush << "\r"; - + #endif + edge_vertices.clear(); get_simplex_vertices( edge, 1, n, binomial_coeff, std::back_inserter(edge_vertices) ); diameters[edge] = dist(edge_vertices[0], edge_vertices[1]); } + #ifdef INDICATE_PROGRESS std::cout << "\033[K"; #endif + #endif // std::cout << "diameters: " << diameters << std::endl; @@ -823,4 +860,4 @@ int main( int argc, char** argv ) { ); -} \ No newline at end of file +} -- cgit v1.2.3