11 #ifndef SPARSE_RIPS_COMPLEX_H_
12 #define SPARSE_RIPS_COMPLEX_H_
14 #include <gudhi/Debug_utils.h>
16 #include <gudhi/choose_n_farthest_points.h>
18 #include <boost/graph/adjacency_list.hpp>
19 #include <boost/range/metafunctions.hpp>
25 namespace rips_complex {
44 template <
typename Filtration_value>
48 typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS,
49 boost::property<vertex_filtration_t, Filtration_value>,
50 boost::property<edge_filtration_t, Filtration_value>>
53 typedef int Vertex_handle;
65 template <
typename RandomAccessPo
intRange,
typename Distance>
68 GUDHI_CHECK(epsilon > 0,
"epsilon must be positive");
69 auto dist_fun = [&](Vertex_handle i, Vertex_handle j) {
return distance(points[i], points[j]); };
71 std::back_inserter(sorted_points), std::back_inserter(params));
72 compute_sparse_graph(dist_fun, epsilon, mini, maxi);
85 template <
typename DistanceMatrix>
88 [&](Vertex_handle i, Vertex_handle j) {
return (i==j) ? 0 : (i<j) ? distance_matrix[j][i] : distance_matrix[i][j]; },
89 epsilon, mini, maxi) {}
101 template <
typename SimplicialComplexForRips>
104 std::invalid_argument(
"Sparse_rips_complex::create_complex - simplicial complex is not empty"));
111 const int n = boost::size(params);
112 std::vector<Filtration_value> lambda(n);
115 lambda[sorted_points[i]] = params[i];
116 double cst = epsilon_ * (1 - epsilon_) / 2;
118 auto filt = complex.filtration(sh);
119 auto mini = filt * cst;
120 for(
auto v : complex.simplex_vertex_range(sh)){
126 complex.expansion_with_blockers(dim_max, block);
131 template <
typename Distance>
133 const auto& points = sorted_points;
134 const int n = boost::size(points);
135 double cst = epsilon * (1 - epsilon) / 2;
137 new (&graph_) Graph(n);
139 typename boost::graph_traits<Graph>::vertex_iterator v_i, v_e;
140 for (std::tie(v_i, v_e) = vertices(graph_); v_i != v_e; ++v_i) {
143 put(vertex_filtration_t(), graph_, v, 0);
149 for (
int i = 0; i < n; ++i) {
150 auto&& pi = points[i];
152 if (li < mini)
break;
153 for (
int j = i + 1; j < n; ++j) {
154 auto&& pj = points[j];
155 auto d = dist(pi, pj);
157 if (lj < mini)
break;
158 GUDHI_CHECK(lj <= li,
"Bad furthest point sorting");
162 if (d * epsilon <= 2 * lj)
164 else if (d * epsilon > li + lj)
167 alpha = (d - lj / epsilon) * 2;
169 if (epsilon < 1 && alpha * cst > lj)
174 add_edge(pi, pj, alpha, graph_);
183 std::vector<Vertex_handle> sorted_points;
185 std::vector<Filtration_value> params;
192 #endif // SPARSE_RIPS_COMPLEX_H_