From bf20ea8c8ff1e272436839f74ae50ad6913e6a95 Mon Sep 17 00:00:00 2001 From: Gard Spreemann Date: Fri, 19 Oct 2018 14:02:22 +0200 Subject: Initial commit of cleaned-up version. --- include/misc.hpp | 68 +++++++++++ include/simplex_tree.hpp | 311 +++++++++++++++++++++++++++++++++++++++++++++++ include/sparse_row.hpp | 48 ++++++++ 3 files changed, 427 insertions(+) create mode 100644 include/misc.hpp create mode 100644 include/simplex_tree.hpp create mode 100644 include/sparse_row.hpp diff --git a/include/misc.hpp b/include/misc.hpp new file mode 100644 index 0000000..614f30f --- /dev/null +++ b/include/misc.hpp @@ -0,0 +1,68 @@ +#pragma once + +#include +#include +#include +#include + +#include + +namespace Gudhi_laplacian +{ + // FIXME: Read from PETSc headers. + static const int32_t PETSC_MAT_FILE_CLASSID = 1211216; + static const int32_t PETSC_VEC_FILE_CLASSID = 1211214; + + bool preexisting_files(const std::string & prefix, const std::string & postfix) + { + glob_t globbuf; + std::string pattern = prefix + std::string(".*.") + postfix; + bool ret = (glob(pattern.c_str(), GLOB_ERR | GLOB_NOSORT, NULL, &globbuf) != GLOB_NOMATCH); + globfree(&globbuf); + return ret; + } + + + template inline void reverse_endianness(T & x) + { + uint8_t * p = reinterpret_cast(&x); + std::reverse(p, p + sizeof(T)); + } + + template inline T read_le(std::istream & s) + { + T result; + s.read(reinterpret_cast(&result), sizeof(T)); + #ifdef BIGENDIAN + reverse_endianness(result); + #endif + return result; + } + + template inline void write_le(std::ostream & s, T value) + { + #ifdef BIGENDIAN + reverse_endianness(value); + #endif + s.write(reinterpret_cast(&value), sizeof(T)); + } + + template inline T read_be(std::istream & s) + { + T result; + s.read(reinterpret_cast(&result), sizeof(T)); + #ifndef BIGENDIAN + reverse_endianness(result); + #endif + return result; + } + + template inline void write_be(std::ostream & s, T value) + { + #ifndef BIGENDIAN + reverse_endianness(value); + #endif + s.write(reinterpret_cast(&value), sizeof(T)); + } +} + diff --git a/include/simplex_tree.hpp b/include/simplex_tree.hpp new file mode 100644 index 0000000..72cae62 --- /dev/null +++ b/include/simplex_tree.hpp @@ -0,0 +1,311 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include + +#include "sparse_row.hpp" +#include "misc.hpp" + +namespace Gudhi_laplacian +{ + typedef Gudhi::Simplex_tree<> ST; + typedef std::vector > Dimensional_stratification; + + // This function is not const in the simplex tree because a bunch + // of boundary/coboundary query functions are not const. + // UB if x is not a face of y. + int face_sign(ST & st, ST::Simplex_handle x, ST::Simplex_handle y) + { + auto x_vers = st.simplex_vertex_range(x); + auto y_vers = st.simplex_vertex_range(y); + + auto it = x_vers.begin(); + auto jt = y_vers.begin(); + + int sign = (st.dimension(y) % 2 == 0) ? 1 : -1; // The vertex ranges run backwards. + while (it != x_vers.end()) + { + if (*it != *jt) + return sign; + + ++it; + ++jt; + sign *= -1; + } + return sign; + } + + Dimensional_stratification make_contiguous_keys(ST & st) + { + Dimensional_stratification stratification; + std::vector keys; + + for (auto s : st.complex_simplex_range()) + { + int dim = st.dimension(s); + + while (stratification.size() <= dim) + stratification.push_back(std::vector()); + while (keys.size() <= dim) + keys.push_back(0); + + stratification[dim].push_back(s); + st.assign_key(s, keys[dim]++); + } + stratification.push_back(std::vector()); + + return stratification; + } + + Sparse_row assemble_cbd_row(ST & st, ST::Simplex_handle s) + { + std::cerr << "WARNING: This function has not been tested much." << std::endl; + Sparse_row 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)); + sign *= -1; + } + + // Note: We don't want to compress the row. + return ret; + } + + Sparse_row assemble_laplacian_row(ST & st, const std::vector > & weights, ST::Simplex_handle s) + { + Sparse_row row; + + int dim = st.dimension(s); + ST::Simplex_key s_key = ST::key(s); + double s_weight = weights[dim][s_key]; + + // Up part. + ST::Cofaces_simplex_range cbds = st.cofaces_simplex_range(s, 1); + for (auto cbd : cbds) + { + int s_cbd_sign = face_sign(st, s, cbd); + ST::Simplex_key cbd_key = ST::key(cbd); + double cbd_weight = weights[dim+1][cbd_key]; + ST::Boundary_simplex_range bds = st.boundary_simplex_range(cbd); + 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))); + bd_cbd_sign *= -1; + } + } + + // Down part. + ST::Boundary_simplex_range bds = st.boundary_simplex_range(s); + int bd_s_sign = (dim % 2 == 0) ? 1 : -1; + for (auto bd : bds) + { + ST::Simplex_key bd_key = ST::key(bd); + double bd_weight = weights[dim-1][bd_key]; + ST::Cofaces_simplex_range cbds = st.cofaces_simplex_range(bd, 1); + for (auto cbd : cbds) + { + 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))); + } + bd_s_sign *= -1; + } + + compress_sparse_row(row); + return row; + } + + int write_cbdrys(ST & st, const Dimensional_stratification & stratification, 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; + + for (int d = from_degree; d < to_degree; ++d) + { + // Coboundary matrix in degree d. + // Shape: (dim C_{d+1}) x (dim C_{d}) + int m = stratification[d+1].size(); + int n = stratification[d].size(); + + 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, m); + write_be(file, n); + write_be(file, (d+2)*m); + + for (int i = 0; i < m; ++i) + write_be(file, d+2); + + // For indices. + for (auto it = stratification[d+1].cbegin(); it != stratification[d+1].cend(); ++it) + { + Sparse_row row = assemble_cbd_row(st, *it); + assert(row.size() == d+2); + for (auto jt = row.cbegin(); jt != row.cend(); ++jt) + { + write_be(file, jt->first); + } + } + + // For values. + for (auto it = stratification[d+1].cbegin(); it != stratification[d+1].cend(); ++it) + { + Sparse_row row = assemble_cbd_row(st, *it); + assert(row.size() == d+2); + for (auto jt = row.cbegin(); jt != row.cend(); ++jt) + { + write_be(file, jt->second); + } + } + + file.close(); + } + + return 0; + } + + int write_laplacians(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; + + for (int d = from_degree; d < to_degree; ++d) + { + // Laplacian matrix in degree d. + // Shape: (dim C_{d}) x (dim C_{d}) + int m = stratification[d].size(); + + 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, m); + write_be(file, m); + + int num_nonzero = 0; + std::vector num_nonzero_row; + num_nonzero_row.reserve(m); + + // First traversal (counting). + for (auto it = stratification[d].cbegin(); it != stratification[d].cend(); ++it) + { + Sparse_row row = assemble_laplacian_row(st, weights, *it); + num_nonzero += row.size(); + num_nonzero_row.push_back(row.size()); + } + + write_be(file, num_nonzero); + for (auto it = num_nonzero_row.cbegin(); it != num_nonzero_row.cend(); ++it) + write_be(file, *it); + + // Second traversal (indices). + for (auto it = stratification[d].cbegin(); it != stratification[d].cend(); ++it) + { + Sparse_row row = assemble_laplacian_row(st, weights, *it); + for (auto jt = row.cbegin(); jt != row.cend(); ++jt) + { + write_be(file, jt->first); + } + } + + // Third traversal (values). + for (auto it = stratification[d].cbegin(); it != stratification[d].cend(); ++it) + { + Sparse_row row = assemble_laplacian_row(st, weights, *it); + for (auto jt = row.cbegin(); jt != row.cend(); ++jt) + { + write_be(file, jt->second); + } + } + + file.close(); + } + + return 0; + } + + int write_weights(const std::vector > & weights, int from_degree, int to_degree, const std::string & prefix) + { + if (preexisting_files(prefix, std::string("petsc"))) + return 1; + + if (weights.size() <= to_degree) + return 2; + + for (int d = from_degree; d < to_degree; ++d) + { + // Weights in degree d. + int m = weights[d].size(); + + 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_VEC_FILE_CLASSID); + write_be(file, m); + + for (auto it = weights[d].cbegin(); it != weights[d].cend(); ++it) + write_be(file, *it); + + file.close(); + } + + return 0; + } + + int write_complex(ST & st, const Dimensional_stratification & stratification, int from_degree, int to_degree, const std::string & prefix) + { + if (preexisting_files(prefix, std::string("txt"))) + return 1; + + if (stratification.size() <= to_degree) + return 2; + + for (int d = from_degree; d < to_degree; ++d) + { + std::ofstream file; + file.open(prefix + std::string(".") + std::to_string(d) + std::string(".txt"), std::ios::out); + + if (!file.is_open()) + return 3; + + for (auto it = stratification[d].begin(); it != stratification[d].end(); ++it) + { + std::vector tmp(std::begin(st.simplex_vertex_range(*it)), std::end(st.simplex_vertex_range(*it))); + std::reverse(std::begin(tmp), std::end(tmp)); + for (auto jt = std::begin(tmp); jt != std::end(tmp) - 1; ++jt) + { + file << std::to_string(*jt) << " "; + } + file << tmp.back() << std::endl; + } + file.close(); + } + + return 0; + } + +} diff --git a/include/sparse_row.hpp b/include/sparse_row.hpp new file mode 100644 index 0000000..9896625 --- /dev/null +++ b/include/sparse_row.hpp @@ -0,0 +1,48 @@ +#pragma once + +#include +#include +#include +#include + +namespace Gudhi_laplacian +{ + template + using Sparse_row = std::vector >; + + template + void compress_sparse_row(Sparse_row & row) + { + static_assert(std::is_arithmetic::value, "Only arithmetic types are supported."); + + Sparse_row tmp(row); + row.clear(); + std::sort(tmp.begin(), tmp.end()); + for (auto it = tmp.cbegin(); it != tmp.cend(); ++it) + { + if (row.empty()) + { + if (it->second != T()) + { + row.push_back(*it); + } + } + else if (it->first == row.back().first) + { + T x = row.back().second + it->second; + if (x == T()) + { + row.pop_back(); + } + else + { + row.back().second = x; + } + } + else if (it->second != T()) + { + row.push_back(*it); + } + } + } +} -- cgit v1.2.3