diff options
author | Marc Glisse <marc.glisse@inria.fr> | 2019-05-10 22:00:54 +0200 |
---|---|---|
committer | Marc Glisse <marc.glisse@inria.fr> | 2019-05-10 22:00:54 +0200 |
commit | fe2f26629481faf316c74dafa7eae892bbc7a556 (patch) | |
tree | 8491fbc96a6fa333ec588bad7352bb508463963e /src/Rips_complex/include/gudhi/Sparse_rips_complex.h | |
parent | d697f28afbcc6be4d5346b6a390b2549d7bf481f (diff) |
Fix blocker formula
I had copied it straight from the paper, without adding the factor 2 we
have everywhere.
Diffstat (limited to 'src/Rips_complex/include/gudhi/Sparse_rips_complex.h')
-rw-r--r-- | src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 13 |
1 files changed, 9 insertions, 4 deletions
diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index a249cd8e..8df6e387 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -126,7 +126,7 @@ class Sparse_rips_complex { // lambda[original_order]=params[sorted_order] for(int i=0;i<n;++i) lambda[sorted_points[i]] = params[i]; - double cst = epsilon_ * (1 - epsilon_); + double cst = epsilon_ * (1 - epsilon_) / 2; auto block = [cst,&complex,&lambda](typename SimplicialComplexForRips::Simplex_handle sh){ auto filt = complex.filtration(sh); auto mini = filt * cst; @@ -156,6 +156,7 @@ class Sparse_rips_complex { void compute_sparse_graph(Distance& dist, double epsilon, Filtration_value mini, Filtration_value maxi) { const auto& points = sorted_points; // convenience alias const int n = boost::size(points); + double cst = epsilon * (1 - epsilon) / 2; graph_.~Graph(); new (&graph_) Graph(n); // for(auto v : vertices(g)) // doesn't work :-( @@ -184,10 +185,14 @@ class Sparse_rips_complex { // The paper has d/2 and d-lj/e to match the Cech, but we use doubles to match the Rips if (d * epsilon <= 2 * lj) alpha = d; - else if (d * epsilon <= li + lj && (epsilon >= 1 || d * epsilon <= lj * (1 + 1 / (1 - epsilon)))) - alpha = (d - lj / epsilon) * 2; - else + else if (d * epsilon > li + lj) continue; + else { + alpha = (d - lj / epsilon) * 2; + // Keep the test exactly the same as in block to avoid inconsistencies + if (epsilon < 1 && alpha * cst > lj) + continue; + } if (alpha <= maxi) add_edge(pi, pj, alpha, graph_); |