From fb018924c2af5a2fa7786e6968b3d638c01e6198 Mon Sep 17 00:00:00 2001 From: Ulrich Bauer Date: Mon, 19 Oct 2015 00:40:42 +0200 Subject: space-efficient persistence computation --- ripser.cpp | 134 +++++++++++++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 122 insertions(+), 12 deletions(-) diff --git a/ripser.cpp b/ripser.cpp index 1740d2c..9129ff3 100644 --- a/ripser.cpp +++ b/ripser.cpp @@ -114,13 +114,13 @@ public: for (int i = 0; i <= dim; ++i) for (int j = i + 1; j <= dim; ++j) { b_diam = std::max(b_diam, dist(vertices[i], vertices[j])); - if (a_diam < b_diam) + 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); // return (a_diam == b_diam) && (((a < b) && !reverse) || ((a > b) && reverse)); } @@ -208,7 +208,29 @@ public: }; +template +int get_pivot(Heap& column) +{ + if( column.empty() ) + return -1; + else { + int max_element = column.top(); + column.pop(); + while( !column.empty() && column.top() == max_element ) { + column.pop(); + if( column.empty() ) + return -1; + else { + max_element = column.top(); + column.pop(); + } + } + if( max_element != -1 ) + column.push( max_element ); + return max_element; + } +} int main() { @@ -234,8 +256,10 @@ int main() { {4,9,7,5,0} }; + int n = dist.size(); + int dim_max = 2; - binomial_coeff_table binomial_coeff(dist.size(), dim_max + 2); + binomial_coeff_table binomial_coeff(n, dim_max + 2); std::cout << dist (0,1) << std::endl; @@ -266,22 +290,23 @@ int main() { int dim = 1; int simplex_index = 2; + double threshold = 7; std::vector vertices; - get_simplex_vertices( simplex_index, dim, dist.size(), binomial_coeff, std::back_inserter(vertices) ); + get_simplex_vertices( simplex_index, dim, n, binomial_coeff, std::back_inserter(vertices) ); std::cout << "coboundary of simplex " << vertices << ":" << std::endl; std::vector coboundary; - get_simplex_coboundary( simplex_index, dim, dist.size(), binomial_coeff, std::back_inserter(coboundary) ); + get_simplex_coboundary( simplex_index, dim, n, binomial_coeff, std::back_inserter(coboundary) ); for (int coboundary_simplex_index: coboundary) { std::vector vertices; - get_simplex_vertices( coboundary_simplex_index, dim + 1, dist.size(), binomial_coeff, std::back_inserter(vertices) ); + get_simplex_vertices( coboundary_simplex_index, dim + 1, n, binomial_coeff, std::back_inserter(vertices) ); std::cout << " " << coboundary_simplex_index << " " << vertices << " (" << comp1.diam(coboundary_simplex_index) << ")" << std::endl; } @@ -310,19 +335,17 @@ int main() { std::cout << "entries: " << csm.entries << std::endl; - std::priority_queue, rips_filtration_comparator > working_boundary(comp1); + std::priority_queue, rips_filtration_comparator > queue(comp1); - for (int e: coboundary) working_boundary.push(e); + for (int e: coboundary) queue.push(e); - std::cout << "pivot of coboundary: " << working_boundary.top() << std::endl; + std::cout << "pivot of coboundary: " << queue.top() << std::endl; std::cout << (comp1(3,6) ? "3<6" : "3>=6") << std::endl; std::cout << (comp1(0,6) ? "0<6" : "0>=6") << std::endl; - std::unordered_map pivot_column_index; - std::vector columns_to_reduce = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 }; std::sort(columns_to_reduce.begin(), columns_to_reduce.end(), comp1); @@ -332,8 +355,95 @@ int main() { for (int index: columns_to_reduce) { std::vector vertices; - get_simplex_vertices( index, dim, dist.size(), binomial_coeff, std::back_inserter(vertices) ); + get_simplex_vertices( index, dim, n, binomial_coeff, std::back_inserter(vertices) ); + std::cout << " " << index << " " << vertices << " (" << comp1.diam(index) << ")" << std::endl; + } + + + + std::unordered_map pivot_column_index; + + for (auto persistence_pair: pivot_column_index) { + + } + + int num_simplices = binomial_coeff(n, dim + 1); + + columns_to_reduce.clear(); + + for (int index = 0; index < num_simplices; ++index) { + if (comp1.diam(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(), comp1); + + std::cout << "sorted columns to reduce: " << columns_to_reduce << std::endl; + + for (int index: columns_to_reduce) { + std::vector vertices; + + get_simplex_vertices( index, dim, n, binomial_coeff, std::back_inserter(vertices) ); std::cout << " " << index << " " << vertices << " (" << comp1.diam(index) << ")" << std::endl; } + + //compressed_sparse_matrix reduction_matrix; + + std::priority_queue, rips_filtration_comparator > working_coboundary(comp2); + + //std::vector reduction_column; + + for (int index: columns_to_reduce) { + //reduction_column.clear(); + + int pivot, column = index; + + do { + + coboundary.clear(); + get_simplex_coboundary( column, dim, n, binomial_coeff, std::back_inserter(coboundary) ); + for (int e: coboundary) if (comp2.diam(e) <= threshold) working_coboundary.push(e); + + 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) { + auto pair = pivot_column_index.find(pivot); + + if (pair == pivot_column_index.end()) { + pivot_column_index.insert(std::make_pair(pivot, index)); + break; + } + + column = pair->second; + } + + /* + coboundary.clear(); + get_simplex_coboundary( birth, dim, n, binomial_coeff, std::back_inserter(coboundary) ); + for (int e: coboundary) if (comp1.diam(e) <= threshold) working_coboundary.push(e); + + //space-efficient variant: drop this part (and the reduction_matrix) + + for (int col = reduction_matrix.cbegin()) { + coboundary.clear(); + get_simplex_coboundary( col, dim, n, binomial_coeff, std::back_inserter(coboundary) ); + for (int e: coboundary) if (comp1.diam(e) <= threshold) working_coboundary.push(e); + } + */ + } while ( pivot != -1 ); + + } + + std::cout << "pairs:" << std::endl << pivot_column_index << std::endl; + + } \ No newline at end of file -- cgit v1.2.3