summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--include/gudhi_laplacian/simplex_tree.hpp127
-rw-r--r--include/gudhi_laplacian/sparse_vector.hpp114
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<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;
}
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 <algorithm>
#include <utility>
#include <type_traits>
+#include <cassert>
+
+#include <iostream>
namespace Gudhi_laplacian
{
template <typename T>
- using Sparse_vector = std::vector<std::pair<unsigned int, T> >;
-
- template <typename T>
- void compress_sparse_vector(Sparse_vector<T> & row)
+ class Sparse_vector
{
- static_assert(std::is_arithmetic<T>::value, "Only arithmetic types are supported.");
-
- Sparse_vector<T> 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<T>::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<std::pair<unsigned int, T> > 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<Sparse_vector<T> > transpose(const std::vector<Sparse_vector<T> > & x, unsigned int n)
+ {
+ unsigned int m = x.size();
+
+ std::vector<Sparse_vector<T> > ret(n);
+ for (unsigned int i = 0; i < m; ++i)
{
- row.push_back(*it);
+ const Sparse_vector<T> & 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<std::pair<unsigned int, T> >::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<unsigned int, T> & x, const std::pair<unsigned int, T> & y)
+ { return x.first < y.first; });
+ if (it == data.cend() || it->first != i)
+ return T();
+ else
+ return it->second;
+ }
+
+ inline typename std::vector<std::pair<unsigned int, T> >::size_type size() const { return data.size(); }
+ inline typename std::vector<std::pair<unsigned int, T> >::size_type nnz() const { return data.size(); }
+
+ using Const_iterator = typename std::vector<std::pair<unsigned int, T> >::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<std::pair<unsigned int, T> > data;
+ };
+
}