From 8b982f5ac72792185d31846df5db3b8ad4d0f2fb Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Thu, 22 Oct 2015 13:49:39 +0200 Subject: cleaned up unused code --- ripser.cpp | 374 ++++++++++--------------------------------------------------- 1 file changed, 62 insertions(+), 312 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 8a5b3ec..c4310c5 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -7,7 +7,6 @@ #include #include #include "prettyprint.hpp" -//#include typedef float value_t; typedef long index_t; @@ -53,15 +52,12 @@ public: } 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); return B[n][k]; } }; - - template OutputIterator get_simplex_vertices( index_t idx, const index_t dim, index_t n, const binomial_coeff_table& binomial_coeff, OutputIterator out ) { @@ -126,17 +122,10 @@ private: public: rips_filtration_comparator( - const DistanceMatrix& _dist, - const index_t _dim, - const binomial_coeff_table& _binomial_coeff -// , -// bool _reverse = true - ): - dist(_dist), dim(_dim), vertices(_dim + 1), - binomial_coeff(_binomial_coeff) -// , -// reverse(_reverse) - {}; + const DistanceMatrix& _dist, + const index_t _dim, + const binomial_coeff_table& _binomial_coeff + ): dist(_dist), dim(_dim), vertices(_dim + 1), binomial_coeff(_binomial_coeff) {}; inline dist_t diameter(const index_t index) const { dist_t diam = 0; @@ -155,7 +144,6 @@ public: assert(b < binomial_coeff(dist.size(), dim + 1)); dist_t a_diam = 0, b_diam = 0; - b_diam = diameter(b); @@ -169,99 +157,8 @@ public: return (a_diam == b_diam) && (a > b); } - }; -template -class rips_filtration_comparator_dim; - -template -class rips_filtration_comparator_dim { -public: - const DistanceMatrix& dist; - -private: - mutable std::vector vertices; - - typedef decltype(dist(0,0)) dist_t; - - bool reverse; - - const binomial_coeff_table& binomial_coeff; - -public: - rips_filtration_comparator_dim( - const DistanceMatrix& _dist, - const binomial_coeff_table& _binomial_coeff - ): - dist(_dist), vertices(3), - binomial_coeff(_binomial_coeff) - {}; - - inline dist_t diameter(const index_t index) const { - dist_t diam = 0; - get_simplex_vertices(index, 2, dist.size(), binomial_coeff, vertices.begin() ); - - diam = std::max(diam, dist(vertices[0], vertices[1])); - diam = std::max(diam, dist(vertices[0], vertices[2])); - diam = std::max(diam, dist(vertices[0], vertices[3])); - diam = std::max(diam, dist(vertices[1], vertices[2])); - diam = std::max(diam, dist(vertices[1], vertices[3])); - diam = std::max(diam, dist(vertices[2], vertices[3])); - - return diam; - } - - 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)); - - dist_t a_diam = 0, b_diam = 0; - - - b_diam = diameter(b); - - get_simplex_vertices(a, 2, dist.size(), binomial_coeff, vertices.begin() ); - if ((a_diam = dist(vertices[0], vertices[1])) >= b_diam) return (a_diam > b_diam) || (a > b); - if ((a_diam = dist(vertices[0], vertices[2])) >= b_diam) return (a_diam > b_diam) || (a > b); - if ((a_diam = dist(vertices[0], vertices[3])) >= b_diam) return (a_diam > b_diam) || (a > b); - if ((a_diam = dist(vertices[1], vertices[2])) >= b_diam) return (a_diam > b_diam) || (a > b); - if ((a_diam = dist(vertices[1], vertices[3])) >= b_diam) return (a_diam > b_diam) || (a > b); - if ((a_diam = dist(vertices[2], vertices[3])) >= b_diam) return (a_diam > b_diam) || (a > b); - - return false; - } - -}; - - -template -void inline get_simplex_coboundary( index_t idx, index_t dim, index_t n, const binomial_coeff_table& binomial_coeff, OutputIterator coboundary ) -{ - - --n; - - index_t modified_idx = idx; - - for( index_t k = dim + 1; k != -1 && n != -1; --k ) { - while( binomial_coeff( n , k ) > idx ) { - *coboundary++ = modified_idx + binomial_coeff( n , k + 1 ); - if (n==0) break; - --n; - } - - idx -= binomial_coeff( n , k ); - - modified_idx -= binomial_coeff( n , k ); - modified_idx += binomial_coeff( n , k + 1 ); - - --n; - } - - return; -} - class simplex_coboundary_enumerator { private: index_t idx, modified_idx, dim, n, k; @@ -270,14 +167,9 @@ private: public: inline simplex_coboundary_enumerator ( - index_t _idx, - index_t _dim, - index_t _n, + index_t _idx, index_t _dim, index_t _n, const binomial_coeff_table& _binomial_coeff): - idx(_idx), modified_idx(_idx), dim(_dim), k(dim + 1), n(_n - 1), binomial_coeff(_binomial_coeff) - { - - } + idx(_idx), modified_idx(_idx), dim(_dim), k(dim + 1), n(_n - 1), binomial_coeff(_binomial_coeff) {} inline bool has_next() { @@ -309,8 +201,7 @@ public: template std::vector get_diameters ( - const DistanceMatrix& dist, - const index_t dim, + const DistanceMatrix& dist, const index_t dim, const std::vector& previous_diameters, const binomial_coeff_table& binomial_coeff ) @@ -339,11 +230,6 @@ std::vector get_diameters ( std::cout << "\033[K"; #endif -// rips_filtration_comparator comp(dist, dim, binomial_coeff); -// for (index_t simplex = 0; simplex < diameters.size(); ++simplex) { -// assert(diameters[simplex] == comp.diameter(simplex)); -// } - return diameters; } @@ -359,7 +245,6 @@ public: typedef value_t dist_t; - const binomial_coeff_table& binomial_coeff; public: @@ -386,7 +271,6 @@ public: return ((a_diam > b_diam) || ((a_diam == b_diam) && (a > b))); } - }; @@ -401,7 +285,6 @@ public: inline size_t size() const { return distances.size(); } - }; @@ -439,7 +322,6 @@ public: inline size_t size() const { return n; } - }; @@ -474,8 +356,6 @@ public: init(); } -// template -// compressed_lower_distance_matrix(const DistanceMatrix& mat) : compressed_lower_distance_matrix(mat.size()) { compressed_lower_distance_matrix(const compressed_upper_distance_matrix& mat) { n = mat.size(); @@ -488,7 +368,6 @@ public: } } - inline value_t operator()(const index_t i, const index_t j) const { if (i > j) return rows[i][j]; @@ -501,44 +380,6 @@ public: inline size_t size() const { return n; } - -}; - -template -class compressed_sparse_matrix { -public: - std::vector bounds; - std::vector entries; - - - inline size_t size() const { - return bounds.size(); - } - - inline typename std::vector::const_iterator cbegin(size_t index) const { - assert(index < size()); - return index == 0 ? entries.cbegin() : entries.cbegin() + bounds[index - 1]; - } - - inline typename std::vector::const_iterator cend(size_t index) const { - assert(index < size()); - return entries.cbegin() + bounds[index]; - } - - template - inline void append(Iterator begin, Iterator end) { - for (Iterator it = begin; it != end; ++it) { - entries.push_back(*it); - } - bounds.push_back(entries.size()); - } - - - template - inline void append(const Collection collection) { - append(collection.cbegin(), collection.cend()); - } - }; template @@ -571,20 +412,6 @@ inline index_t get_pivot(Heap& column) return max_element; } -template -inline std::vector move_to_column_vector(Heap& column) -{ - std::vector temp_col; - index_t pivot = pop_pivot( column ); - while( pivot != -1 ) { - temp_col.push_back( pivot ); - pivot = pop_pivot( column ); - } - return temp_col; -} - - - void print_help_and_exit() { std::cerr << "Usage: " << "ripser " << "[options] input_filename output_filename" << std::endl; @@ -608,85 +435,73 @@ void compute_pairs( const binomial_coeff_table& binomial_coeff ) { - std::cout << "persistence intervals in dim " << dim << ":" << std::endl; + 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]; + + #ifdef ASSEMBLE_REDUCTION_COLUMN + std::priority_queue, decltype(comp) > reduction_column(comp); + #endif + + 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; + + std::vector coboundary; + + do { - for (index_t i = 0; i < columns_to_reduce.size(); ++i) { - index_t index = columns_to_reduce[i]; - #ifdef ASSEMBLE_REDUCTION_COLUMN - std::priority_queue, decltype(comp) > reduction_column(comp); - #endif - - 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"; + reduction_column.push( column ); #endif - - index_t pivot, column = index; - - std::vector coboundary; - do { + simplex_coboundary_enumerator coboundaries(column, dim, n, binomial_coeff); + while (coboundaries.has_next()) { + index_t coface = coboundaries.next(); + if (comp.diameter(coface) <= threshold) + working_coboundary.push(coface); + } + + pivot = get_pivot(working_coboundary); - #ifdef ASSEMBLE_REDUCTION_COLUMN - reduction_column.push( column ); - #endif - - simplex_coboundary_enumerator coboundaries(column, dim, n, binomial_coeff); - while (coboundaries.has_next()) { - index_t coface = coboundaries.next(); - if (comp.diameter(coface) <= threshold) - working_coboundary.push(coface); - } - - pivot = get_pivot(working_coboundary); + if (pivot != -1) { + auto pair = pivot_column_index.find(pivot); - - //add boundary column at index birth to working_coboundary - - //since the boundary column is not stored explicitly, - //add the boundary of the column at index birth in the reduction matrix instead - - //space-efficient variant: add just the boundary column at index birth instead - //this avoids having to store the reduction matrix - - if (pivot != -1) { - auto pair = pivot_column_index.find(pivot); - - 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; - } + if (pair == pivot_column_index.end()) { + pivot_column_index.insert(std::make_pair(pivot, index)); - column = pair->second; - } + #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 - } 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; + break; + } + + column = pair->second; } - #endif + } 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; } - - std::cout << "\033[K"; - -} + #endif + } + std::cout << "\033[K"; +} int main( int argc, char** argv ) { @@ -740,7 +555,6 @@ int main( int argc, char** argv ) { exit(-1); } - int64_t file_type; input_stream.read( (char*)&file_type, sizeof( int64_t ) ); if( file_type != 7 ) { @@ -751,18 +565,12 @@ int main( int argc, char** argv ) { int64_t n; input_stream.read( (char*)&n, sizeof( int64_t ) ); - - //std::vector distances = - distance_matrix dist; dist.distances = std::vector>(n, std::vector(n)); for( int i = 0; i < n; ++i ) { input_stream.read( (char*)&dist.distances[i][0], n * sizeof(int64_t) ); } - - //std::cout << dist.distances << std::endl; - #endif @@ -773,8 +581,6 @@ int main( int argc, char** argv ) { while(std::getline(input_stream, value_string, ',')) distances.push_back(std::stod(value_string)); -// std::cout << "distances: " << distances << std::endl; - compressed_upper_distance_matrix dist_upper(std::move(distances)); index_t n = dist_upper.size(); @@ -782,44 +588,15 @@ int main( int argc, char** argv ) { #endif std::cout << "distance matrix with " << n << " points" << std::endl; - -// std::vector> distance_matrix_full(n, std::vector(n)); -// for (index_t i = 0; i < n; ++i) -// for (index_t j = 0; j < n; ++j) -// distance_matrix_full[i][j] = dist_upper(i, j); -// std::cout << distance_matrix_full << std::endl; - - - + compressed_lower_distance_matrix dist(dist_upper); std::cout << "distance matrix transformed to lower triangular form" << std::endl; -// for (index_t i = 0; i < n; ++i) -// for (index_t j = 0; j < n; ++j) -// distance_matrix_full[i][j] = dist(i, j); -// std::cout << distance_matrix_full << std::endl; - - -// std::cout << "distances (lower): " << dist.distances << std::endl; - -// return 0; - assert(dim_max < n - 2); binomial_coeff_table binomial_coeff(n, dim_max + 2); -// std::vector coboundary; -// get_simplex_coboundary( 1, 2, n, binomial_coeff, std::back_inserter(coboundary) ); -// std::cout << coboundary < columns_to_reduce; @@ -833,30 +610,7 @@ int main( int argc, char** argv ) { std::vector& diameters = dist.distances; - -// std::cout << "precomputing 1-simplex diameters" << std::endl; -// -// std::vector diameters( binomial_coeff( n, 2 ) ); -// -// 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; - for (index_t dim = 0; dim < dim_max; ++dim) { @@ -879,7 +633,6 @@ int main( int argc, char** argv ) { binomial_coeff ); - index_t num_simplices = binomial_coeff(n, dim + 2); columns_to_reduce.clear(); @@ -893,8 +646,6 @@ int main( int argc, char** argv ) { std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), comp); - - #ifdef PRECOMPUTE_DIAMETERS previous_diameters = std::move(diameters); @@ -906,7 +657,6 @@ int main( int argc, char** argv ) { std::cout << "precomputing " << dim + 2 << "-simplex diameters" << std::endl; diameters = get_diameters( dist, dim + 2, previous_diameters, binomial_coeff ); #endif - } index_t dim = dim_max; -- cgit v1.2.3