diff options
Diffstat (limited to 'include/gudhi_laplacian/simplex_tree.hpp')
-rw-r--r-- | include/gudhi_laplacian/simplex_tree.hpp | 127 |
1 files changed, 122 insertions, 5 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<int> assemble_cbd_row(ST & st, ST::Simplex_handle s) { - std::cerr << "WARNING: This function has not been tested much." << std::endl; Sparse_vector<int> 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<std::vector<double> > & 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<Sparse_vector<int> > cbd_same; + std::vector<Sparse_vector<int> > cbd_same_t; + std::vector<Sparse_vector<int> > 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<Sparse_vector<int> >(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<int>::transpose(cbd_same, n); + + std::vector<Sparse_vector<double> > lap; + unsigned int nnz = 0; + + for (unsigned int i = 0; i < n; ++i) + { + lap.push_back(Sparse_vector<double>()); + + // 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<int32_t>(file, PETSC_MAT_FILE_CLASSID); + write_be<int32_t>(file, n); + write_be<int32_t>(file, n); + write_be<int32_t>(file, nnz); + for (unsigned int i = 0; i < n; ++i) + write_be<int32_t>(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<int32_t>(file, it->first); + } + } + for (unsigned int i = 0; i < n; ++i) + { + for (auto it = lap[i].cbegin(); it != lap[i].cend(); ++it) + { + write_be<double>(file, it->second); + } + } + + file.close(); + } + return 0; + } Sparse_vector<double> assemble_laplacian_row(ST & st, const std::vector<std::vector<double> > & 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; } |