From 05be96e7cf650913075cfcbb114287520afcbca6 Mon Sep 17 00:00:00 2001 From: Gard Spreemann Date: Wed, 7 Nov 2018 14:52:41 +0100 Subject: Rewrite of sparse vector + inefficient mm-based Laplacian. --- include/gudhi_laplacian/simplex_tree.hpp | 127 ++++++++++++++++++++++++++++-- include/gudhi_laplacian/sparse_vector.hpp | 114 +++++++++++++++++++++------ 2 files changed, 212 insertions(+), 29 deletions(-) diff --git a/include/gudhi_laplacian/simplex_tree.hpp b/include/gudhi_laplacian/simplex_tree.hpp index dafcd06..6d5aff5 100644 --- a/include/gudhi_laplacian/simplex_tree.hpp +++ b/include/gudhi_laplacian/simplex_tree.hpp @@ -64,20 +64,137 @@ namespace Gudhi_laplacian Sparse_vector assemble_cbd_row(ST & st, ST::Simplex_handle s) { - std::cerr << "WARNING: This function has not been tested much." << std::endl; Sparse_vector ret; // This is a specialization of face_sign. int sign = (st.dimension(s) % 2 == 0) ? 1 : -1; for (auto bd : st.boundary_simplex_range(s)) { - ret.push_back(std::make_pair(ST::key(bd), sign)); + ret.add(ST::key(bd), sign); sign *= -1; } // Note: We don't want to compress the row. return ret; } + + int write_laplacians_mm(ST & st, const Dimensional_stratification & stratification, const std::vector > & weights, int from_degree, int to_degree, const std::string & prefix) + { + if (preexisting_files(prefix, std::string("petsc"))) + return 1; + + if (stratification.size() <= to_degree) + return 2; + + std::vector > cbd_same; + std::vector > cbd_same_t; + std::vector > cbd_below; + + for (int d = from_degree; d < to_degree; ++d) + { + unsigned int m = d+1 < stratification.size() ? stratification[d+1].size() : 0; + unsigned int n = stratification[d].size(); + unsigned int o = d > 0 ? stratification[d-1].size() : 0; + + if (d == 0) + { + cbd_below = std::vector >(n); // An nx0 matrix to avoid special-casing further down. + } + else + { + cbd_below = cbd_same; + } + assert(cbd_below.size() == n); + + cbd_same.clear(); + for (unsigned int i = 0; i < m; ++i) + { + cbd_same.push_back(assemble_cbd_row(st, stratification[d+1][i])); + } + cbd_same_t = Sparse_vector::transpose(cbd_same, n); + + std::vector > lap; + unsigned int nnz = 0; + + for (unsigned int i = 0; i < n; ++i) + { + lap.push_back(Sparse_vector()); + + // Up part. + for (unsigned int j = 0; j < n; ++j) + { + auto it = cbd_same_t[i].cbegin(); + auto jt = cbd_same_t[j].cbegin(); + while (it != cbd_same_t[i].cend() && jt != cbd_same_t[j].cend()) + { + if (it->first < jt->first) + ++it; + else if (jt->first < it->first) + ++jt; + else + { + double foo = (weights[d+1][jt->first]/weights[d][i])*(it->second * jt->second); + lap[i].add(j, (weights[d+1][jt->first]/weights[d][i])*(it->second * jt->second)); + ++it; + ++jt; + } + } + } + + // Down part. + for (unsigned int j = 0; j < n; ++j) + { + auto it = cbd_below[i].cbegin(); + auto jt = cbd_below[j].cbegin(); + while (it != cbd_below[i].cend() && jt != cbd_below[j].cend()) + { + if (it->first < jt->first) + ++it; + else if (jt->first < it->first) + ++jt; + else + { + lap[i].add(j, (weights[d][j]/weights[d-1][it->first])*(it->second * jt->second)); + ++it; + ++jt; + } + } + } + + lap[i].compress(); + nnz += lap[i].nnz(); + } + + std::ofstream file; + file.open(prefix + std::string(".") + std::to_string(d) + std::string(".petsc"), std::ios::out | std::ios::binary); + if (!file.is_open()) + return 3; + + write_be(file, PETSC_MAT_FILE_CLASSID); + write_be(file, n); + write_be(file, n); + write_be(file, nnz); + for (unsigned int i = 0; i < n; ++i) + write_be(file, lap[i].nnz()); + for (unsigned int i = 0; i < n; ++i) + { + for (auto it = lap[i].cbegin(); it != lap[i].cend(); ++it) + { + write_be(file, it->first); + } + } + for (unsigned int i = 0; i < n; ++i) + { + for (auto it = lap[i].cbegin(); it != lap[i].cend(); ++it) + { + write_be(file, it->second); + } + } + + file.close(); + } + return 0; + } Sparse_vector assemble_laplacian_row(ST & st, const std::vector > & weights, ST::Simplex_handle s) { @@ -98,7 +215,7 @@ namespace Gudhi_laplacian int bd_cbd_sign = ((dim+1) % 2 == 0) ? 1 : -1; for (auto bd : bds) { - row.push_back(std::make_pair(ST::key(bd), s_cbd_sign * bd_cbd_sign * (cbd_weight/s_weight))); + row.add(ST::key(bd), s_cbd_sign * bd_cbd_sign * (cbd_weight/s_weight)); bd_cbd_sign *= -1; } } @@ -116,12 +233,12 @@ namespace Gudhi_laplacian int bd_cbd_sign = face_sign(st, bd, cbd); ST::Simplex_key cbd_key = ST::key(cbd); double cbd_weight = weights[dim][cbd_key]; - row.push_back(std::make_pair(cbd_key, bd_s_sign * bd_cbd_sign * (cbd_weight/bd_weight))); + row.add(cbd_key, bd_s_sign * bd_cbd_sign * (cbd_weight/bd_weight)); } bd_s_sign *= -1; } - compress_sparse_vector(row); + row.compress(); return row; } diff --git a/include/gudhi_laplacian/sparse_vector.hpp b/include/gudhi_laplacian/sparse_vector.hpp index 58dd40c..4fff58a 100644 --- a/include/gudhi_laplacian/sparse_vector.hpp +++ b/include/gudhi_laplacian/sparse_vector.hpp @@ -4,45 +4,111 @@ #include #include #include +#include + +#include namespace Gudhi_laplacian { template - using Sparse_vector = std::vector >; - - template - void compress_sparse_vector(Sparse_vector & row) + class Sparse_vector { - static_assert(std::is_arithmetic::value, "Only arithmetic types are supported."); - - Sparse_vector tmp(row); - row.clear(); - std::sort(tmp.begin(), tmp.end()); - for (auto it = tmp.cbegin(); it != tmp.cend(); ++it) + public: + Sparse_vector() : compressed(true), data() + { + static_assert(std::is_arithmetic::value, "Only arithmetic types are supported."); + } + + // Default destructor. + // Default copy-constructor. + // Default assignment operator. + // Default moves. + + void compress() { - if (row.empty()) + if (compressed) + return; + + std::vector > tmp(data); + data.clear(); + std::sort(tmp.begin(), tmp.end()); + for (auto it = tmp.cbegin(); it != tmp.cend(); ++it) { - if (it->second != T()) + if (data.empty()) { - row.push_back(*it); + if (it->second != T()) + { + data.push_back(*it); + } } - } - else if (it->first == row.back().first) - { - T x = row.back().second + it->second; - if (x == T()) + else if (it->first == data.back().first) { - row.pop_back(); + T x = data.back().second + it->second; + if (x == T()) + { + data.pop_back(); + } + else + { + data.back().second = x; + } } - else + else if (it->second != T()) { - row.back().second = x; + data.push_back(*it); } } - else if (it->second != T()) + compressed = true; + } + + inline void add(unsigned int i, T x) + { + data.push_back(std::make_pair(i, x)); + compressed = false; + } + + static std::vector > transpose(const std::vector > & x, unsigned int n) + { + unsigned int m = x.size(); + + std::vector > ret(n); + for (unsigned int i = 0; i < m; ++i) { - row.push_back(*it); + const Sparse_vector & a = x[i]; + for (auto it = a.cbegin(); it != a.cend(); ++it) + { + assert(it->first < n); + ret[it->first].data.push_back(std::make_pair(i, it->second)); // Avoid add to keep compressed = true. + } } + return ret; } - } + + inline const T & operator[](typename std::vector >::size_type i) + { + std::cerr << "UNTESTED FUNCTION" << std::endl; + assert(compressed); + auto it = std::lower_bound(data.cbegin(), data.cend(), std::make_pair(i, T()), + [](const std::pair & x, const std::pair & y) + { return x.first < y.first; }); + if (it == data.cend() || it->first != i) + return T(); + else + return it->second; + } + + inline typename std::vector >::size_type size() const { return data.size(); } + inline typename std::vector >::size_type nnz() const { return data.size(); } + + using Const_iterator = typename std::vector >::const_iterator; + inline Const_iterator cbegin() const { return data.cbegin(); } + inline Const_iterator cend() const { return data.cend(); } + inline Const_iterator begin() const { return data.cbegin(); } + inline Const_iterator end() const { return data.cend(); } + + private: + bool compressed; + std::vector > data; + }; + } -- cgit v1.2.3