From b159c58b09a7385aedf645f5086232267837b459 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Wed, 21 Oct 2015 10:26:03 +0200 Subject: more options for precomputed diameters --- ripser.cpp | 463 +++++++++++++++++++++++++++++++++++-------------------------- 1 file changed, 265 insertions(+), 198 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 4ab08d9..0a62675 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -10,7 +10,10 @@ //#include -//#define PRECOMPUTE_DIAMETERS +#define PRECOMPUTE_DIAMETERS + +#define PRECOMPUTE_DIAMETERS_IN_TOP_DIMENSION + //#define ASSEMBLE_REDUCTION_COLUMN @@ -43,7 +46,8 @@ public: long operator()(long n, long k) const { // std::cout << "B(" << n << "," << k << ")\n"; - + assert(n <= n_max); + assert(k <= k_max); return B[n][k]; } }; @@ -110,7 +114,6 @@ public: for (long i = 0; i <= dim; ++i) for (long j = i + 1; j <= dim; ++j) { - auto d = dist(vertices[i], vertices[j]); diam = std::max(diam, dist(vertices[i], vertices[j])); } return diam; @@ -231,12 +234,12 @@ public: binomial_coeff(_binomial_coeff) {}; - double diameter(const long a) { + double diameter(const long a) const { assert(a < diameters.size()); return diameters[a]; } - bool operator()(const long a, const long b) + bool operator()(const long a, const long b) const { assert(a < diameters.size()); assert(b < diameters.size()); @@ -395,6 +398,216 @@ void print_help_and_exit() exit(-1); } +template +void compute_pairs( + std::vector& columns_to_reduce, + const ComparatorCofaces& comp, const Comparator& comp_prev, + long dim, long dim_max, long n, + double threshold, + const binomial_coeff_table& binomial_coeff +) { + std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), comp_prev); + + + std::unordered_map pivot_column_index; + + std::cout << "persistence intervals in dim " << dim << ":" << std::endl; + +// std::cout << "reduce columns in dim " << dim << ": " << columns_to_reduce << std::endl; +// std::cout << " rows in dim " << dim + 1 << ": " << columns_to_reduce << std::endl; + +// std::vector rows(binomial_coeff(n, dim + 2)); +// for (long simplex_index = 0; simplex_index < binomial_coeff(n, dim + 2); ++simplex_index) { +// rows[simplex_index] = simplex_index; +// } +// std::sort(rows.begin(), rows.end(), comp); +// +// for (long simplex_index: rows) { +// std::vector vertices; +// +// get_simplex_vertices( simplex_index, dim + 1, n, binomial_coeff, std::back_inserter(vertices) ); +// std::cout << " " << simplex_index << " " << vertices << " (" << comp.diameter(simplex_index) << ")" << std::endl; +// +// } + + for (long i = 0; i < columns_to_reduce.size(); ++i) { + long 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); + +// std::cout << "reduce column " << index << std::setw(0) +// << " (" << i + 1 << "/" << columns_to_reduce.size() << ")\r"; + + std::cout << "\033[K" << "reducing column " << i + 1 << "/" << columns_to_reduce.size() + << " (diameter " << comp_prev.diameter(index) << ")" + << std::flush << "\r"; + + + long pivot, column = index; + + std::vector coboundary; + + + do { + + #ifdef ASSEMBLE_REDUCTION_COLUMN + reduction_column.push( column ); + #endif + + coboundary.clear(); + get_simplex_coboundary( column, dim, n, binomial_coeff, std::back_inserter(coboundary) ); + +// std::vector sorted_coboundary = coboundary; std::sort(sorted_coboundary.begin(), sorted_coboundary.end(), comp); +// std::cout << "add " << sorted_coboundary << " to working col" << std::endl; +// for (long coboundary_simplex_index: coboundary) { +// std::vector vertices; +// +// get_simplex_vertices( coboundary_simplex_index, dim + 1, n, binomial_coeff, std::back_inserter(vertices) ); +// std::cout << " " << coboundary_simplex_index << " " << vertices << " (" << comp.diameter(coboundary_simplex_index) << ")" << std::endl; +// } + + for (long e: coboundary) if (comp.diameter(e) <= threshold) working_coboundary.push(e); // std::cout << "push " << e << std::endl;} + +// std::cout << "=" << std::endl; +// auto working_coboundary_copy = working_coboundary; +// while (!working_coboundary_copy.empty()) { +// std::cout << " " << working_coboundary_copy.top() << std::endl; +// working_coboundary_copy.pop(); +// } +// std::cout << "=" << std::endl; + + + pivot = get_pivot(working_coboundary); + + + //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) { +// std::cout << "pivot: " << pivot << std::endl; + auto pair = pivot_column_index.find(pivot); + + if (pair == pivot_column_index.end()) { + pivot_column_index.insert(std::make_pair(pivot, index)); + + double birth = comp_prev.diameter(index), death = comp.diameter(pivot); + if (birth != death) + std::cout << "\033[K" << " [" << birth << "," << death << ")" << std::endl << std::flush; + + break; + } + + column = pair->second; + } + + /* + coboundary.clear(); + get_simplex_coboundary( birth, dim, n, binomial_coeff, std::back_inserter(coboundary) ); + for (long e: coboundary) if (comp2.diameter(e) <= threshold) working_coboundary.push(e); + + //space-efficient variant: drop this part (and the reduction_matrix) + + for (long col = reduction_matrix.cbegin()) { + coboundary.clear(); + get_simplex_coboundary( col, dim, n, binomial_coeff, std::back_inserter(coboundary) ); + for (long e: coboundary) if (comp2.diameter(e) <= threshold) working_coboundary.push(e); + } + */ + } while ( pivot != -1 ); + + + if ( pivot == -1 ) { + double birth = comp_prev.diameter(index); + std::cout << "\033[K" << " [" << birth << ", )" << std::endl << std::flush; + } + + +// std::cout << "size of working column heap: " << working_coboundary.size() << std::endl; +// +// #ifdef ASSEMBLE_REDUCTION_COLUMN +// std::vector cochain = move_to_column_vector( reduction_column ); +//// std::cout << "reduction cochain: " << cochain << std::endl; +// +// if ( pivot != -1 ) { +// std::vector coboundary = move_to_column_vector( working_coboundary ); +//// std::cout << "reduced coboundary: " << coboundary << std::endl; +// } +// +// std::cout << "fill-in: " << cochain.size()-1 << std::endl << std::endl; +// #endif + + + } + + +// std::cout << "dimension " << dim << " pairs:" << std::endl; +//// std::cout << pivot_column_index << std::endl; +// +// +// +// +// for (std::pair pair: pivot_column_index) { +// double birth = comp_prev.diameter(pair.second), death = comp.diameter(pair.first); +//// std::cout << vertices_of_simplex(pair.second, dim, n, binomial_coeff) << "," << +//// vertices_of_simplex(pair.first, dim + 1, n, binomial_coeff) << std::endl; +// if (birth != death) +// std::cout << " [" << birth << "," << death << ")" << std::endl; +// } + + if (dim == dim_max) + return; + + + long num_simplices = binomial_coeff(n, dim + 2); + + columns_to_reduce.clear(); + +// std::cout << "columns to reduce in dim " << dim + 1 << " (" << num_simplices << " total)" << std::endl; + + for (long index = 0; index < num_simplices; ++index) { + +// if (comp.diameter(index) > threshold) { +// std::cout << " " << vertices_of_simplex(index, dim + 1, n, binomial_coeff) << ": " << comp.diameter(index) << " above threshold" << std::endl; +// } else if (pivot_column_index.find(index) != pivot_column_index.end()) { +// std::cout << " " << vertices_of_simplex(index, dim + 1, n, binomial_coeff) << " appears in pair" << std::endl; +// } + + if (comp.diameter(index) <= threshold && pivot_column_index.find(index) == pivot_column_index.end()) { + columns_to_reduce.push_back(index); + } + } + + + + + +// std::cout << "sorted " << dim + 1 << "-columns to reduce: " << columns_to_reduce << std::endl; +// +// for (long index: columns_to_reduce) { +// std::vector vertices; +// +// get_simplex_vertices( index, dim, n, binomial_coeff, std::back_inserter(vertices) ); +// std::cout << " " << index << " " << vertices << " (" << comp.diameter(index) << ")" << std::endl; +// } +// +// std::cout << std::endl; + + + +} + + + + int main( int argc, char** argv ) { if( argc < 3 ) print_help_and_exit(); @@ -402,7 +615,7 @@ int main( int argc, char** argv ) { std::string input_filename = argv[ argc - 2 ]; std::string output_filename = argv[ argc - 1 ]; - long dim_max = 2; + long dim_max = 1; double threshold = std::numeric_limits::max(); for( long idx = 1; idx < argc - 2; idx++ ) { @@ -514,7 +727,7 @@ int main( int argc, char** argv ) { // n = dist.size(); - assert(dim_max < n - 1); + assert(dim_max < n - 2); binomial_coeff_table binomial_coeff(n, dim_max + 2); @@ -635,7 +848,6 @@ int main( int argc, char** argv ) { std::vector columns_to_reduce; - std::vector coboundary; for (long index = n - 1; index >= 0; --index) { columns_to_reduce.push_back(index); @@ -654,218 +866,73 @@ int main( int argc, char** argv ) { for (long edge = 0; edge < diameters.size(); ++edge) { // diameters[edge] = comp1.diameter(edge); + std::cout << "\033[Kstoring diameter of simplex " << edge + 1 << "/" << diameters.size() << std::flush << "\r"; + 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]); } + std::cout << "\033[K"; #endif + + + for (long dim = 0; dim < dim_max; ++dim) { - - - //compressed_sparse_matrix reduction_matrix; - - + #ifdef PRECOMPUTE_DIAMETERS rips_filtration_diameter_comparator comp(diameters, dim + 1, binomial_coeff); rips_filtration_diameter_comparator comp_prev(previous_diameters, dim, binomial_coeff); #else rips_filtration_comparator comp(dist, dim + 1, binomial_coeff); rips_filtration_comparator comp_prev(dist, dim, binomial_coeff); - #endif - - - std::unordered_map pivot_column_index; - - std::cout << "persistence intervals in dim " << dim << ":" << std::endl; - -// std::cout << "reduce columns in dim " << dim << ": " << columns_to_reduce << std::endl; -// std::cout << " rows in dim " << dim + 1 << ": " << columns_to_reduce << std::endl; - -// std::vector rows(binomial_coeff(n, dim + 2)); -// for (long simplex_index = 0; simplex_index < binomial_coeff(n, dim + 2); ++simplex_index) { -// rows[simplex_index] = simplex_index; -// } -// std::sort(rows.begin(), rows.end(), comp); -// -// for (long simplex_index: rows) { -// std::vector vertices; -// -// get_simplex_vertices( simplex_index, dim + 1, n, binomial_coeff, std::back_inserter(vertices) ); -// std::cout << " " << simplex_index << " " << vertices << " (" << comp.diameter(simplex_index) << ")" << std::endl; -// -// } - - for (long i = 0; i < columns_to_reduce.size(); ++i) { - long 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); - -// std::cout << "reduce column " << index << std::setw(0) -// << " (" << i + 1 << "/" << columns_to_reduce.size() << ")\r"; - - std::cout << "\033[K" << "reducing column " << i + 1 << "/" << columns_to_reduce.size() - << " (diameter " << comp_prev.diameter(index) << ")" - << std::flush << "\r"; - - - long pivot, column = index; - - do { - - #ifdef ASSEMBLE_REDUCTION_COLUMN - reduction_column.push( column ); - #endif - - coboundary.clear(); - get_simplex_coboundary( column, dim, n, binomial_coeff, std::back_inserter(coboundary) ); - - std::vector sorted_coboundary = coboundary; std::sort(sorted_coboundary.begin(), sorted_coboundary.end(), comp); -// std::cout << "add " << sorted_coboundary << " to working col" << std::endl; -// for (long coboundary_simplex_index: coboundary) { -// std::vector vertices; -// -// get_simplex_vertices( coboundary_simplex_index, dim + 1, n, binomial_coeff, std::back_inserter(vertices) ); -// std::cout << " " << coboundary_simplex_index << " " << vertices << " (" << comp.diameter(coboundary_simplex_index) << ")" << std::endl; -// } - - for (long e: coboundary) if (comp.diameter(e) <= threshold) working_coboundary.push(e); // std::cout << "push " << e << std::endl;} - -// std::cout << "=" << std::endl; -// auto working_coboundary_copy = working_coboundary; -// while (!working_coboundary_copy.empty()) { -// std::cout << " " << working_coboundary_copy.top() << std::endl; -// working_coboundary_copy.pop(); -// } -// std::cout << "=" << std::endl; - - - pivot = get_pivot(working_coboundary); - - - //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) { -// std::cout << "pivot: " << pivot << std::endl; - auto pair = pivot_column_index.find(pivot); - - if (pair == pivot_column_index.end()) { - pivot_column_index.insert(std::make_pair(pivot, index)); - - double birth = comp_prev.diameter(index), death = comp.diameter(pivot); - if (birth != death) - std::cout << "\033[K" << " [" << birth << "," << death << ")" << std::endl << std::flush; - - break; - } - - column = pair->second; - } - - /* - coboundary.clear(); - get_simplex_coboundary( birth, dim, n, binomial_coeff, std::back_inserter(coboundary) ); - for (long e: coboundary) if (comp2.diameter(e) <= threshold) working_coboundary.push(e); - - //space-efficient variant: drop this part (and the reduction_matrix) - - for (long col = reduction_matrix.cbegin()) { - coboundary.clear(); - get_simplex_coboundary( col, dim, n, binomial_coeff, std::back_inserter(coboundary) ); - for (long e: coboundary) if (comp2.diameter(e) <= threshold) working_coboundary.push(e); - } - */ - } while ( pivot != -1 ); - - -// std::cout << "size of working column heap: " << working_coboundary.size() << std::endl; -// -// #ifdef ASSEMBLE_REDUCTION_COLUMN -// std::vector cochain = move_to_column_vector( reduction_column ); -//// std::cout << "reduction cochain: " << cochain << std::endl; -// -// if ( pivot != -1 ) { -// std::vector coboundary = move_to_column_vector( working_coboundary ); -//// std::cout << "reduced coboundary: " << coboundary << std::endl; -// } -// -// std::cout << "fill-in: " << cochain.size()-1 << std::endl << std::endl; -// #endif - - - } - - -// std::cout << "dimension " << dim << " pairs:" << std::endl; -//// std::cout << pivot_column_index << std::endl; -// -// -// -// -// for (std::pair pair: pivot_column_index) { -// double birth = comp_prev.diameter(pair.second), death = comp.diameter(pair.first); -//// std::cout << vertices_of_simplex(pair.second, dim, n, binomial_coeff) << "," << -//// vertices_of_simplex(pair.first, dim + 1, n, binomial_coeff) << std::endl; -// if (birth != death) -// std::cout << " [" << birth << "," << death << ")" << std::endl; -// } - - if (dim == dim_max - 1) - break; - - - long num_simplices = binomial_coeff(n, dim + 2); - - columns_to_reduce.clear(); - -// std::cout << "columns to reduce in dim " << dim + 1 << " (" << num_simplices << " total)" << std::endl; - - for (long index = 0; index < num_simplices; ++index) { - -// if (comp.diameter(index) > threshold) { -// std::cout << " " << vertices_of_simplex(index, dim + 1, n, binomial_coeff) << ": " << comp.diameter(index) << " above threshold" << std::endl; -// } else if (pivot_column_index.find(index) != pivot_column_index.end()) { -// std::cout << " " << vertices_of_simplex(index, dim + 1, n, binomial_coeff) << " appears in pair" << std::endl; -// } - - if (comp.diameter(index) <= threshold && pivot_column_index.find(index) == pivot_column_index.end()) { - columns_to_reduce.push_back(index); - } - } - std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), comp); + #endif -// std::cout << "sorted " << dim + 1 << "-columns to reduce: " << columns_to_reduce << std::endl; -// -// for (long index: columns_to_reduce) { -// std::vector vertices; -// -// get_simplex_vertices( index, dim, n, binomial_coeff, std::back_inserter(vertices) ); -// std::cout << " " << index << " " << vertices << " (" << comp.diameter(index) << ")" << std::endl; -// } -// -// std::cout << std::endl; - + compute_pairs( + columns_to_reduce, + comp, comp_prev, + dim, dim_max, n, + threshold, + binomial_coeff + ); std::cout << "\033[K"; #ifdef PRECOMPUTE_DIAMETERS previous_diameters = std::move(diameters); - std::cout << "precomputing " << dim + 1 << "-simplex diameters" << std::endl; + + #ifndef PRECOMPUTE_DIAMETERS_IN_TOP_DIMENSION + if (dim == dim_max - 1) + break; + #endif + + std::cout << "precomputing " << dim + 2 << "-simplex diameters" << std::endl; diameters = get_diameters( dist, dim + 2, previous_diameters, binomial_coeff ); #endif - - } + long dim = dim_max; + + #ifdef PRECOMPUTE_DIAMETERS + rips_filtration_diameter_comparator comp_prev(previous_diameters, dim, binomial_coeff); + #else + rips_filtration_comparator comp_prev(dist, dim, binomial_coeff); + #endif + + #ifdef PRECOMPUTE_DIAMETERS_IN_TOP_DIMENSION + rips_filtration_diameter_comparator comp(diameters, dim + 1, binomial_coeff); + #else + rips_filtration_comparator comp(dist, dim + 1, binomial_coeff); + #endif + + compute_pairs( + columns_to_reduce, + comp, comp_prev, + dim, dim_max, n, + threshold, + binomial_coeff + ); + + } \ No newline at end of file -- cgit v1.2.3