diff options
Diffstat (limited to 'src/Rips_complex/include/gudhi/Sparse_rips_complex.h')
-rw-r--r-- | src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 79 |
1 files changed, 60 insertions, 19 deletions
diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 4dcc08ed..8df6e387 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -47,7 +47,9 @@ namespace rips_complex { * * \details * This class is used to construct a sparse \f$(1+O(\epsilon))\f$-approximation of `Rips_complex`, i.e. a filtered - * simplicial complex that is multiplicatively \f$(1+O(\epsilon))\f$-interleaved with the Rips filtration. + * simplicial complex that is multiplicatively + * \f$(1+O(\epsilon))\f$-interleaved with the Rips filtration. More precisely, + * this is a \f$(1,\frac{1}{1-\epsilon}\f$-interleaving. * * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex. */ @@ -55,7 +57,7 @@ template <typename Filtration_value> class Sparse_rips_complex { private: // TODO(MG): use a different graph where we know we can safely insert in parallel. - typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS, + typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS, boost::property<vertex_filtration_t, Filtration_value>, boost::property<edge_filtration_t, Filtration_value>> Graph; @@ -68,32 +70,36 @@ class Sparse_rips_complex { * @param[in] points Range of points. * @param[in] distance Distance function that returns a `Filtration_value` from 2 given points. * @param[in] epsilon Approximation parameter. epsilon must be positive. + * @param[in] mini Minimal filtration value. Ignore anything below this scale. This is a less efficient version of `Gudhi::subsampling::sparsify_point_set()`. + * @param[in] maxi Maximal filtration value. Ignore anything above this scale. * */ template <typename RandomAccessPointRange, typename Distance> - Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) { + Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon, Filtration_value mini=-std::numeric_limits<Filtration_value>::infinity(), Filtration_value maxi=std::numeric_limits<Filtration_value>::infinity()) + : epsilon_(epsilon) { GUDHI_CHECK(epsilon > 0, "epsilon must be positive"); - std::vector<Vertex_handle> sorted_points; - std::vector<Filtration_value> params; auto dist_fun = [&](Vertex_handle i, Vertex_handle j) { return distance(points[i], points[j]); }; Ker<decltype(dist_fun)> kernel(dist_fun); subsampling::choose_n_farthest_points(kernel, boost::irange<Vertex_handle>(0, boost::size(points)), -1, -1, std::back_inserter(sorted_points), std::back_inserter(params)); - compute_sparse_graph(sorted_points, params, dist_fun, epsilon); + compute_sparse_graph(dist_fun, epsilon, mini, maxi); } /** \brief Sparse_rips_complex constructor from a distance matrix. * * @param[in] distance_matrix Range of range of distances. * `distance_matrix[i][j]` returns the distance between points \f$i\f$ and - * \f$j\f$ as long as \f$ 0 \leqslant i < j \leqslant + * \f$j\f$ as long as \f$ 0 \leqslant j < i \leqslant * distance\_matrix.size().\f$ * @param[in] epsilon Approximation parameter. epsilon must be positive. + * @param[in] mini Minimal filtration value. Ignore anything below this scale. This is a less efficient version of `Gudhi::subsampling::sparsify_point_set()`. + * @param[in] maxi Maximal filtration value. Ignore anything above this scale. */ template <typename DistanceMatrix> - Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon) + Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon, Filtration_value mini=-std::numeric_limits<Filtration_value>::infinity(), Filtration_value maxi=std::numeric_limits<Filtration_value>::infinity()) : Sparse_rips_complex(boost::irange<Vertex_handle>(0, boost::size(distance_matrix)), - [&](Vertex_handle i, Vertex_handle j) { return distance_matrix[j][i]; }, epsilon) {} + [&](Vertex_handle i, Vertex_handle j) { return (i==j) ? 0 : (i<j) ? distance_matrix[j][i] : distance_matrix[i][j]; }, + epsilon, mini, maxi) {} /** \brief Fills the simplicial complex with the sparse Rips graph and * expands it with all the cliques, stopping at a given maximal dimension. @@ -111,7 +117,26 @@ class Sparse_rips_complex { std::invalid_argument("Sparse_rips_complex::create_complex - simplicial complex is not empty")); complex.insert_graph(graph_); - complex.expansion(dim_max); + if(epsilon_ >= 1) { + complex.expansion(dim_max); + return; + } + const int n = boost::size(params); + std::vector<Filtration_value> lambda(n); + // lambda[original_order]=params[sorted_order] + for(int i=0;i<n;++i) + lambda[sorted_points[i]] = params[i]; + 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; + for(auto v : complex.simplex_vertex_range(sh)){ + if(lambda[v] < mini) + return true; // v died before this simplex could be born + } + return false; + }; + complex.expansion_with_blockers(dim_max, block); } private: @@ -127,9 +152,11 @@ class Sparse_rips_complex { }; // PointRange must be random access. - template <typename PointRange, typename ParamRange, typename Distance> - void compute_sparse_graph(const PointRange& points, const ParamRange& params, Distance& dist, double epsilon) { + template <typename Distance> + 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 :-( @@ -143,29 +170,43 @@ class Sparse_rips_complex { // TODO(MG): // - make it parallel // - only test near-enough neighbors - for (int i = 0; i < n; ++i) + for (int i = 0; i < n; ++i) { + auto&& pi = points[i]; + auto li = params[i]; + if (li < mini) break; for (int j = i + 1; j < n; ++j) { - auto&& pi = points[i]; auto&& pj = points[j]; auto d = dist(pi, pj); - auto li = params[i]; auto lj = params[j]; + if (lj < mini) break; GUDHI_CHECK(lj <= li, "Bad furthest point sorting"); Filtration_value alpha; // 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; + } - add_edge(pi, pj, alpha, graph_); + if (alpha <= maxi) + add_edge(pi, pj, alpha, graph_); } + } } Graph graph_; + double epsilon_; + // Because of the arbitrary split between constructor and create_complex + // sorted_points[sorted_order]=original_order + std::vector<Vertex_handle> sorted_points; + // params[sorted_order]=distance to previous points + std::vector<Filtration_value> params; }; } // namespace rips_complex |