summaryrefslogtreecommitdiff
path: root/include
diff options
context:
space:
mode:
Diffstat (limited to 'include')
-rw-r--r--include/gudhi_laplacian/simplex_tree.hpp46
1 files changed, 17 insertions, 29 deletions
diff --git a/include/gudhi_laplacian/simplex_tree.hpp b/include/gudhi_laplacian/simplex_tree.hpp
index 6d5aff5..179362f 100644
--- a/include/gudhi_laplacian/simplex_tree.hpp
+++ b/include/gudhi_laplacian/simplex_tree.hpp
@@ -1,5 +1,10 @@
#pragma once
+
+/*
+WARNING: DO NOT USE FOR WEIGHTED STUFF YET! ALL WEIGHTS MUST BE 1!
+*/
+
#include <vector>
#include <algorithm>
#include <cassert>
@@ -89,6 +94,7 @@ namespace Gudhi_laplacian
std::vector<Sparse_vector<int> > cbd_same;
std::vector<Sparse_vector<int> > cbd_same_t;
std::vector<Sparse_vector<int> > cbd_below;
+ std::vector<Sparse_vector<int> > cbd_below_t;
for (int d = from_degree; d < to_degree; ++d)
{
@@ -99,12 +105,15 @@ namespace Gudhi_laplacian
if (d == 0)
{
cbd_below = std::vector<Sparse_vector<int> >(n); // An nx0 matrix to avoid special-casing further down.
+ cbd_below_t = Sparse_vector<int>::transpose(cbd_below, 0);
}
else
{
cbd_below = cbd_same;
+ cbd_below_t = cbd_same_t;
}
assert(cbd_below.size() == n);
+ assert(cbd_below_t.size() == o);
cbd_same.clear();
for (unsigned int i = 0; i < m; ++i)
@@ -112,6 +121,8 @@ namespace Gudhi_laplacian
cbd_same.push_back(assemble_cbd_row(st, stratification[d+1][i]));
}
cbd_same_t = Sparse_vector<int>::transpose(cbd_same, n);
+ assert(cbd_same.size() == m);
+ assert(cbd_same_t.size() == n);
std::vector<Sparse_vector<double> > lap;
unsigned int nnz = 0;
@@ -121,43 +132,20 @@ namespace Gudhi_laplacian
lap.push_back(Sparse_vector<double>());
// Up part.
- for (unsigned int j = 0; j < n; ++j)
+ for (auto it = cbd_same_t[i].cbegin(); it != cbd_same_t[i].cend(); ++it)
{
- 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())
+ for (auto jt = cbd_same[it->first].cbegin(); jt != cbd_same[it->first].cend(); ++jt)
{
- 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;
- }
+ lap[i].add(jt->first, (weights[d+1][it->first]/weights[d][jt->first]) * (it->second * jt->second));
}
}
// Down part.
- for (unsigned int j = 0; j < n; ++j)
+ for (auto it = cbd_below[i].cbegin(); it != cbd_below[i].cend(); ++it)
{
- auto it = cbd_below[i].cbegin();
- auto jt = cbd_below[j].cbegin();
- while (it != cbd_below[i].cend() && jt != cbd_below[j].cend())
+ for (auto jt = cbd_below_t[it->first].cbegin(); jt != cbd_below_t[it->first].cend(); ++jt)
{
- 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].add(jt->first, (weights[d][jt->first]/weights[d-1][it->first]) * (it->second * jt->second));
}
}