23 #ifndef TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_ 24 #define TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_ 26 #include <gudhi/Tangential_complex/config.h> 27 #include <gudhi/Tangential_complex/utilities.h> 28 #include <gudhi/Debug_utils.h> 29 #include <gudhi/console_color.h> 31 #include <CGAL/iterator.h> 34 #include <boost/graph/graph_traits.hpp> 35 #include <boost/graph/adjacency_list.hpp> 36 #include <boost/graph/connected_components.hpp> 37 #include <boost/container/flat_set.hpp> 47 namespace tangential_complex {
50 class Simplicial_complex {
52 typedef boost::container::flat_set<std::size_t> Simplex;
53 typedef std::set<Simplex> Simplex_set;
60 const Simplex &s,
bool perform_checks =
true) {
62 unsigned int num_pts =
static_cast<int> (s.size());
63 std::vector<Complex::iterator> to_erase;
64 bool check_higher_dim_simpl =
true;
65 for (Complex::iterator it_simplex = m_complex.begin(),
66 it_simplex_end = m_complex.end();
67 it_simplex != it_simplex_end;
70 if (check_higher_dim_simpl
71 && it_simplex->size() > num_pts
72 && std::includes(it_simplex->begin(), it_simplex->end(),
73 s.begin(), s.end())) {
78 if (it_simplex->size() < num_pts
79 && std::includes(s.begin(), s.end(),
80 it_simplex->begin(), it_simplex->end())) {
81 to_erase.push_back(it_simplex);
83 check_higher_dim_simpl =
false;
86 for (std::vector<Complex::iterator>::const_iterator it = to_erase.begin();
87 it != to_erase.end(); ++it) {
91 return m_complex.insert(s).second;
94 const Simplex_set &simplex_range()
const {
99 return m_complex.empty();
106 template <
typename Test,
typename Output_it>
107 void get_simplices_matching_test(Test test, Output_it out) {
108 for (Complex::const_iterator it_simplex = m_complex.begin(),
109 it_simplex_end = m_complex.end();
110 it_simplex != it_simplex_end;
112 if (test(*it_simplex))
113 *out++ = *it_simplex;
120 void collapse(
int max_simplex_dim,
bool quiet =
false) {
123 std::cerr <<
"Collapsing... ";
126 int k = max_simplex_dim - 1;
128 typedef Complex::iterator Simplex_iterator;
129 typedef std::vector<Simplex_iterator> Simplex_iterator_list;
130 typedef std::map<Simplex, Simplex_iterator_list> Cofaces_map;
132 std::size_t num_collapsed_maximal_simplices = 0;
134 num_collapsed_maximal_simplices = 0;
137 Cofaces_map cofaces_map;
138 for (Complex::const_iterator it_simplex = m_complex.begin(),
139 it_simplex_end = m_complex.end();
140 it_simplex != it_simplex_end;
142 if (static_cast<int> (it_simplex->size()) > k + 1) {
143 std::vector<Simplex> k_faces;
145 combinations(*it_simplex, k + 1, std::back_inserter(k_faces));
146 for (
const auto &comb : k_faces)
147 cofaces_map[comb].push_back(it_simplex);
158 for (Cofaces_map::const_iterator it_map_elt = cofaces_map.begin(),
159 it_map_end = cofaces_map.end();
160 it_map_elt != it_map_end;
162 if (it_map_elt->second.size() == 1) {
163 std::vector<Simplex> k_faces;
164 const Simplex_iterator_list::value_type &it_Cf =
165 *it_map_elt->second.begin();
166 GUDHI_CHECK(it_Cf->size() == max_simplex_dim + 1,
167 std::logic_error(
"Wrong dimension"));
169 combinations(*it_Cf, k + 1, std::back_inserter(k_faces));
170 for (
const auto &f2 : k_faces) {
172 if (f2 != it_map_elt->first) {
173 Cofaces_map::iterator it_comb_in_map = cofaces_map.find(f2);
174 if (it_comb_in_map->second.size() == 1) {
175 it_comb_in_map->second.clear();
176 m_complex.insert(f2);
178 Simplex_iterator_list::iterator it = std::find(it_comb_in_map->second.begin(),
179 it_comb_in_map->second.end(),
181 GUDHI_CHECK(it != it_comb_in_map->second.end(),
182 std::logic_error(
"Error: it == it_comb_in_map->second.end()"));
183 it_comb_in_map->second.erase(it);
187 m_complex.erase(it_Cf);
188 ++num_collapsed_maximal_simplices;
192 }
while (num_collapsed_maximal_simplices > 0);
196 collapse(max_simplex_dim - 1,
true);
200 std::cerr <<
"done.\n";
204 void display_stats()
const {
205 std::cerr << yellow <<
"Complex stats:\n" << white;
207 if (m_complex.empty()) {
208 std::cerr <<
" * No simplices.\n";
211 std::map<int, std::size_t> simplex_stats;
213 for (Complex::const_iterator it_simplex = m_complex.begin(),
214 it_simplex_end = m_complex.end();
215 it_simplex != it_simplex_end;
217 ++simplex_stats[
static_cast<int> (it_simplex->size()) - 1];
220 for (std::map<int, std::size_t>::const_iterator it_map = simplex_stats.begin();
221 it_map != simplex_stats.end(); ++it_map) {
222 std::cerr <<
" * " << it_map->first <<
"-simplices: " 223 << it_map->second <<
"\n";
229 bool is_pure_pseudomanifold__do_not_check_if_stars_are_connected(
int simplex_dim,
230 bool allow_borders =
false,
231 bool exit_at_the_first_problem =
false,
232 int verbose_level = 0,
233 std::size_t *p_num_wrong_dim_simplices = NULL,
234 std::size_t *p_num_wrong_number_of_cofaces = NULL)
const {
235 typedef Simplex K_1_face;
236 typedef std::map<K_1_face, std::size_t> Cofaces_map;
238 std::size_t num_wrong_dim_simplices = 0;
239 std::size_t num_wrong_number_of_cofaces = 0;
245 Cofaces_map cofaces_map;
246 for (Complex::const_iterator it_simplex = m_complex.begin(),
247 it_simplex_end = m_complex.end();
248 it_simplex != it_simplex_end;
250 if (static_cast<int> (it_simplex->size()) != simplex_dim + 1) {
251 if (verbose_level >= 2)
252 std::cerr <<
"Found a simplex with dim = " 253 << it_simplex->size() - 1 <<
"\n";
254 ++num_wrong_dim_simplices;
256 std::vector<K_1_face> k_1_faces;
259 *it_simplex, simplex_dim, std::back_inserter(k_1_faces));
260 for (
const auto &k_1_face : k_1_faces) {
261 ++cofaces_map[k_1_face];
266 for (Cofaces_map::const_iterator it_map_elt = cofaces_map.begin(),
267 it_map_end = cofaces_map.end();
268 it_map_elt != it_map_end;
270 if (it_map_elt->second != 2
271 && (!allow_borders || it_map_elt->second != 1)) {
272 if (verbose_level >= 2)
273 std::cerr <<
"Found a k-1-face with " 274 << it_map_elt->second <<
" cofaces\n";
276 if (exit_at_the_first_problem)
279 ++num_wrong_number_of_cofaces;
283 bool ret = num_wrong_dim_simplices == 0 && num_wrong_number_of_cofaces == 0;
285 if (verbose_level >= 1) {
286 std::cerr <<
"Pure pseudo-manifold: ";
288 std::cerr << green <<
"YES" << white <<
"\n";
290 std::cerr << red <<
"NO" << white <<
"\n" 291 <<
" * Number of wrong dimension simplices: " 292 << num_wrong_dim_simplices <<
"\n" 293 <<
" * Number of wrong number of cofaces: " 294 << num_wrong_number_of_cofaces <<
"\n";
298 if (p_num_wrong_dim_simplices)
299 *p_num_wrong_dim_simplices = num_wrong_dim_simplices;
300 if (p_num_wrong_number_of_cofaces)
301 *p_num_wrong_number_of_cofaces = num_wrong_number_of_cofaces;
307 std::size_t num_K_simplices()
const {
308 Simplex_set k_simplices;
310 for (Complex::const_iterator it_simplex = m_complex.begin(),
311 it_simplex_end = m_complex.end();
312 it_simplex != it_simplex_end;
314 if (it_simplex->size() == K + 1) {
315 k_simplices.insert(*it_simplex);
316 }
else if (it_simplex->size() > K + 1) {
319 *it_simplex, K + 1, std::inserter(k_simplices, k_simplices.begin()));
323 return k_simplices.size();
326 std::ptrdiff_t euler_characteristic(
bool verbose =
false)
const {
328 std::cerr <<
"\nComputing Euler characteristic of the complex...\n";
330 std::size_t num_vertices = num_K_simplices<0>();
331 std::size_t num_edges = num_K_simplices<1>();
332 std::size_t num_triangles = num_K_simplices<2>();
335 (std::ptrdiff_t) num_vertices
336 - (std::ptrdiff_t) num_edges
337 + (std::ptrdiff_t) num_triangles;
340 std::cerr <<
"Euler characteristic: V - E + F = " 341 << num_vertices <<
" - " << num_edges <<
" + " << num_triangles <<
" = " 351 bool is_pure_pseudomanifold(
353 std::size_t num_vertices,
354 bool allow_borders =
false,
355 bool exit_at_the_first_problem =
false,
356 int verbose_level = 0,
357 std::size_t *p_num_wrong_dim_simplices = NULL,
358 std::size_t *p_num_wrong_number_of_cofaces = NULL,
359 std::size_t *p_num_unconnected_stars = NULL,
360 Simplex_set *p_wrong_dim_simplices = NULL,
361 Simplex_set *p_wrong_number_of_cofaces_simplices = NULL,
362 Simplex_set *p_unconnected_stars_simplices = NULL)
const {
364 if (simplex_dim == 1) {
365 if (p_num_unconnected_stars)
366 *p_num_unconnected_stars = 0;
367 return is_pure_pseudomanifold__do_not_check_if_stars_are_connected(simplex_dim,
369 exit_at_the_first_problem,
371 p_num_wrong_dim_simplices,
372 p_num_wrong_number_of_cofaces);
376 typedef std::vector<std::vector<Complex::const_iterator> > Stars;
377 std::size_t num_wrong_dim_simplices = 0;
378 std::size_t num_wrong_number_of_cofaces = 0;
379 std::size_t num_unconnected_stars = 0;
383 stars.resize(num_vertices);
384 for (Complex::const_iterator it_simplex = m_complex.begin(),
385 it_simplex_end = m_complex.end();
386 it_simplex != it_simplex_end;
388 if (static_cast<int> (it_simplex->size()) != simplex_dim + 1) {
389 if (verbose_level >= 2)
390 std::cerr <<
"Found a simplex with dim = " 391 << it_simplex->size() - 1 <<
"\n";
392 ++num_wrong_dim_simplices;
393 if (p_wrong_dim_simplices)
394 p_wrong_dim_simplices->insert(*it_simplex);
396 for (Simplex::const_iterator it_point_idx = it_simplex->begin();
397 it_point_idx != it_simplex->end();
399 stars[*it_point_idx].push_back(it_simplex);
408 std::size_t center_vertex_index = 0;
409 for (Stars::const_iterator it_star = stars.begin();
410 it_star != stars.end();
411 ++it_star, ++center_vertex_index) {
412 typedef std::map<Simplex, std::vector<std::size_t> >
413 Dm1_faces_to_adj_D_faces;
414 Dm1_faces_to_adj_D_faces dm1_faces_to_adj_d_faces;
416 for (std::size_t i_dsimpl = 0; i_dsimpl < it_star->size(); ++i_dsimpl) {
417 Simplex dm1_simpl_of_link = *((*it_star)[i_dsimpl]);
418 dm1_simpl_of_link.erase(center_vertex_index);
420 std::vector<std::size_t> dm1_simpl_of_link_vec(
421 dm1_simpl_of_link.begin(), dm1_simpl_of_link.end());
423 CGAL::Combination_enumerator<int> dm2_simplices(
424 simplex_dim - 1, 0, simplex_dim);
425 for (; !dm2_simplices.finished(); ++dm2_simplices) {
427 for (
int j = 0; j < simplex_dim - 1; ++j)
428 dm2_simpl.insert(dm1_simpl_of_link_vec[dm2_simplices[j]]);
429 dm1_faces_to_adj_d_faces[dm2_simpl].push_back(i_dsimpl);
434 std::vector<Graph_vertex> d_faces_descriptors;
435 d_faces_descriptors.resize(it_star->size());
436 for (std::size_t j = 0; j < it_star->size(); ++j)
437 d_faces_descriptors[j] = boost::add_vertex(adj_graph);
439 Dm1_faces_to_adj_D_faces::const_iterator dm1_to_d_it =
440 dm1_faces_to_adj_d_faces.begin();
441 Dm1_faces_to_adj_D_faces::const_iterator dm1_to_d_it_end =
442 dm1_faces_to_adj_d_faces.end();
443 for (std::size_t i_km1_face = 0;
444 dm1_to_d_it != dm1_to_d_it_end;
445 ++dm1_to_d_it, ++i_km1_face) {
446 Graph_vertex km1_gv = boost::add_vertex(adj_graph);
448 for (std::vector<std::size_t>::const_iterator kface_it =
449 dm1_to_d_it->second.begin();
450 kface_it != dm1_to_d_it->second.end();
452 boost::add_edge(km1_gv, *kface_it, adj_graph);
455 if (dm1_to_d_it->second.size() != 2
456 && (!allow_borders || dm1_to_d_it->second.size() != 1)) {
457 ++num_wrong_number_of_cofaces;
458 if (p_wrong_number_of_cofaces_simplices) {
459 for (
auto idx : dm1_to_d_it->second)
460 p_wrong_number_of_cofaces_simplices->insert(*((*it_star)[idx]));
466 bool is_connected =
true;
467 if (boost::num_vertices(adj_graph) > 0) {
468 std::vector<int> components(boost::num_vertices(adj_graph));
470 (boost::connected_components(adj_graph, &components[0]) == 1);
474 if (verbose_level >= 2)
475 std::cerr <<
"Error: star #" << center_vertex_index
476 <<
" is not connected\n";
477 ++num_unconnected_stars;
478 if (p_unconnected_stars_simplices) {
479 for (std::vector<Complex::const_iterator>::const_iterator
480 it_simpl = it_star->begin(),
481 it_simpl_end = it_star->end();
482 it_simpl != it_simpl_end;
484 p_unconnected_stars_simplices->insert(**it_simpl);
491 num_wrong_number_of_cofaces /= simplex_dim;
494 num_wrong_dim_simplices == 0
495 && num_wrong_number_of_cofaces == 0
496 && num_unconnected_stars == 0;
498 if (verbose_level >= 1) {
499 std::cerr <<
"Pure pseudo-manifold: ";
501 std::cerr << green <<
"YES" << white <<
"\n";
503 std::cerr << red <<
"NO" << white <<
"\n" 504 <<
" * Number of wrong dimension simplices: " 505 << num_wrong_dim_simplices <<
"\n" 506 <<
" * Number of wrong number of cofaces: " 507 << num_wrong_number_of_cofaces <<
"\n" 508 <<
" * Number of not-connected stars: " 509 << num_unconnected_stars <<
"\n";
513 if (p_num_wrong_dim_simplices)
514 *p_num_wrong_dim_simplices = num_wrong_dim_simplices;
515 if (p_num_wrong_number_of_cofaces)
516 *p_num_wrong_number_of_cofaces = num_wrong_number_of_cofaces;
517 if (p_num_unconnected_stars)
518 *p_num_unconnected_stars = num_unconnected_stars;
524 typedef Simplex_set Complex;
527 typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Adj_graph;
529 typedef boost::graph_traits<Adj_graph>::vertex_descriptor Graph_vertex;
530 typedef boost::graph_traits<Adj_graph>::edge_descriptor Graph_edge;
539 #endif // TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_ Definition: SimplicialComplexForAlpha.h:26