Simplicial_complex.h
1 /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2  * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3  * Author(s): Clement Jamin
4  *
5  * Copyright (C) 2016 Inria
6  *
7  * Modification(s):
8  * - YYYY/MM Author: Description of the modification
9  */
10 
11 #ifndef TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_
12 #define TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_
13 
14 #include <gudhi/Tangential_complex/config.h>
15 #include <gudhi/Tangential_complex/utilities.h>
16 #include <gudhi/Debug_utils.h>
17 #include <gudhi/console_color.h>
18 
19 #include <CGAL/iterator.h>
20 
21 // For is_pure_pseudomanifold
22 #include <boost/graph/graph_traits.hpp>
23 #include <boost/graph/adjacency_list.hpp>
24 #include <boost/graph/connected_components.hpp>
25 #include <boost/container/flat_set.hpp>
26 
27 #include <algorithm>
28 #include <string>
29 #include <fstream>
30 #include <map> // for map<>
31 #include <vector> // for vector<>
32 #include <set> // for set<>
33 
34 namespace Gudhi {
35 namespace tangential_complex {
36 namespace internal {
37 
38 class Simplicial_complex {
39  public:
40  typedef boost::container::flat_set<std::size_t> Simplex;
41  typedef std::set<Simplex> Simplex_set;
42 
43  // If perform_checks = true, the function:
44  // - won't insert the simplex if it is already in a higher dim simplex
45  // - will erase any lower-dim simplices that are faces of the new simplex
46  // Returns true if the simplex was added
47  bool add_simplex(
48  const Simplex &s, bool perform_checks = true) {
49  if (perform_checks) {
50  unsigned int num_pts = static_cast<int> (s.size());
51  std::vector<Complex::iterator> to_erase;
52  bool check_higher_dim_simpl = true;
53  for (Complex::iterator it_simplex = m_complex.begin(),
54  it_simplex_end = m_complex.end();
55  it_simplex != it_simplex_end;
56  ++it_simplex) {
57  // Check if the simplex is not already in a higher dim simplex
58  if (check_higher_dim_simpl
59  && it_simplex->size() > num_pts
60  && std::includes(it_simplex->begin(), it_simplex->end(),
61  s.begin(), s.end())) {
62  // No need to insert it, then
63  return false;
64  }
65  // Check if the simplex includes some lower-dim simplices
66  if (it_simplex->size() < num_pts
67  && std::includes(s.begin(), s.end(),
68  it_simplex->begin(), it_simplex->end())) {
69  to_erase.push_back(it_simplex);
70  // We don't need to check higher-sim simplices any more
71  check_higher_dim_simpl = false;
72  }
73  }
74  for (std::vector<Complex::iterator>::const_iterator it = to_erase.begin();
75  it != to_erase.end(); ++it) {
76  m_complex.erase(*it);
77  }
78  }
79  return m_complex.insert(s).second;
80  }
81 
82  const Simplex_set &simplex_range() const {
83  return m_complex;
84  }
85 
86  bool empty() {
87  return m_complex.empty();
88  }
89 
90  void clear() {
91  m_complex.clear();
92  }
93 
94  template <typename Test, typename Output_it>
95  void get_simplices_matching_test(Test test, Output_it out) {
96  for (Complex::const_iterator it_simplex = m_complex.begin(),
97  it_simplex_end = m_complex.end();
98  it_simplex != it_simplex_end;
99  ++it_simplex) {
100  if (test(*it_simplex))
101  *out++ = *it_simplex;
102  }
103  }
104 
105  // When a simplex S has only one co-face C, we can remove S and C
106  // without changing the topology
107 
108  void collapse(int max_simplex_dim, bool quiet = false) {
109 #ifdef DEBUG_TRACES
110  if (!quiet)
111  std::cerr << "Collapsing... ";
112 #endif
113  // We note k = max_simplex_dim - 1
114  int k = max_simplex_dim - 1;
115 
116  typedef Complex::iterator Simplex_iterator;
117  typedef std::vector<Simplex_iterator> Simplex_iterator_list;
118  typedef std::map<Simplex, Simplex_iterator_list> Cofaces_map;
119 
120  std::size_t num_collapsed_maximal_simplices = 0;
121  do {
122  num_collapsed_maximal_simplices = 0;
123  // Create a map associating each non-maximal k-faces to the list of its
124  // maximal cofaces
125  Cofaces_map cofaces_map;
126  for (Complex::const_iterator it_simplex = m_complex.begin(),
127  it_simplex_end = m_complex.end();
128  it_simplex != it_simplex_end;
129  ++it_simplex) {
130  if (static_cast<int> (it_simplex->size()) > k + 1) {
131  std::vector<Simplex> k_faces;
132  // Get the k-faces composing the simplex
133  combinations(*it_simplex, k + 1, std::back_inserter(k_faces));
134  for (const auto &comb : k_faces)
135  cofaces_map[comb].push_back(it_simplex);
136  }
137  }
138 
139  // For each non-maximal k-face F, if F has only one maximal coface Cf:
140  // - Look for the other k-faces F2, F3... of Cf in the map and:
141  // * if the list contains only Cf, clear the list (we don't remove the
142  // list since it creates troubles with the iterators) and add the F2,
143  // F3... to the complex
144  // * otherwise, remove Cf from the associated list
145  // - Remove Cf from the complex
146  for (Cofaces_map::const_iterator it_map_elt = cofaces_map.begin(),
147  it_map_end = cofaces_map.end();
148  it_map_elt != it_map_end;
149  ++it_map_elt) {
150  if (it_map_elt->second.size() == 1) {
151  std::vector<Simplex> k_faces;
152  const Simplex_iterator_list::value_type &it_Cf =
153  *it_map_elt->second.begin();
154  GUDHI_CHECK(it_Cf->size() == max_simplex_dim + 1,
155  std::logic_error("Wrong dimension"));
156  // Get the k-faces composing the simplex
157  combinations(*it_Cf, k + 1, std::back_inserter(k_faces));
158  for (const auto &f2 : k_faces) {
159  // Skip F
160  if (f2 != it_map_elt->first) {
161  Cofaces_map::iterator it_comb_in_map = cofaces_map.find(f2);
162  if (it_comb_in_map->second.size() == 1) {
163  it_comb_in_map->second.clear();
164  m_complex.insert(f2);
165  } else { // it_comb_in_map->second.size() > 1
166  Simplex_iterator_list::iterator it = std::find(it_comb_in_map->second.begin(),
167  it_comb_in_map->second.end(),
168  it_Cf);
169  GUDHI_CHECK(it != it_comb_in_map->second.end(),
170  std::logic_error("Error: it == it_comb_in_map->second.end()"));
171  it_comb_in_map->second.erase(it);
172  }
173  }
174  }
175  m_complex.erase(it_Cf);
176  ++num_collapsed_maximal_simplices;
177  }
178  }
179  // Repeat until no maximal simplex got removed
180  } while (num_collapsed_maximal_simplices > 0);
181 
182  // Collapse the lower dimension simplices
183  if (k > 0)
184  collapse(max_simplex_dim - 1, true);
185 
186 #ifdef DEBUG_TRACES
187  if (!quiet)
188  std::cerr << "done.\n";
189 #endif
190  }
191 
192  void display_stats() const {
193  std::cerr << yellow << "Complex stats:\n" << white;
194 
195  if (m_complex.empty()) {
196  std::cerr << " * No simplices.\n";
197  } else {
198  // Number of simplex for each dimension
199  std::map<int, std::size_t> simplex_stats;
200 
201  for (Complex::const_iterator it_simplex = m_complex.begin(),
202  it_simplex_end = m_complex.end();
203  it_simplex != it_simplex_end;
204  ++it_simplex) {
205  ++simplex_stats[static_cast<int> (it_simplex->size()) - 1];
206  }
207 
208  for (std::map<int, std::size_t>::const_iterator it_map = simplex_stats.begin();
209  it_map != simplex_stats.end(); ++it_map) {
210  std::cerr << " * " << it_map->first << "-simplices: "
211  << it_map->second << "\n";
212  }
213  }
214  }
215 
216  // verbose_level = 0, 1 or 2
217  bool is_pure_pseudomanifold__do_not_check_if_stars_are_connected(int simplex_dim,
218  bool allow_borders = false,
219  bool exit_at_the_first_problem = false,
220  int verbose_level = 0,
221  std::size_t *p_num_wrong_dim_simplices = NULL,
222  std::size_t *p_num_wrong_number_of_cofaces = NULL) const {
223  typedef Simplex K_1_face;
224  typedef std::map<K_1_face, std::size_t> Cofaces_map;
225 
226  std::size_t num_wrong_dim_simplices = 0;
227  std::size_t num_wrong_number_of_cofaces = 0;
228 
229  // Counts the number of cofaces of each K_1_face
230 
231  // Create a map associating each non-maximal k-faces to the list of its
232  // maximal cofaces
233  Cofaces_map cofaces_map;
234  for (Complex::const_iterator it_simplex = m_complex.begin(),
235  it_simplex_end = m_complex.end();
236  it_simplex != it_simplex_end;
237  ++it_simplex) {
238  if (static_cast<int> (it_simplex->size()) != simplex_dim + 1) {
239  if (verbose_level >= 2)
240  std::cerr << "Found a simplex with dim = "
241  << it_simplex->size() - 1 << "\n";
242  ++num_wrong_dim_simplices;
243  } else {
244  std::vector<K_1_face> k_1_faces;
245  // Get the facets composing the simplex
246  combinations(
247  *it_simplex, simplex_dim, std::back_inserter(k_1_faces));
248  for (const auto &k_1_face : k_1_faces) {
249  ++cofaces_map[k_1_face];
250  }
251  }
252  }
253 
254  for (Cofaces_map::const_iterator it_map_elt = cofaces_map.begin(),
255  it_map_end = cofaces_map.end();
256  it_map_elt != it_map_end;
257  ++it_map_elt) {
258  if (it_map_elt->second != 2
259  && (!allow_borders || it_map_elt->second != 1)) {
260  if (verbose_level >= 2)
261  std::cerr << "Found a k-1-face with "
262  << it_map_elt->second << " cofaces\n";
263 
264  if (exit_at_the_first_problem)
265  return false;
266  else
267  ++num_wrong_number_of_cofaces;
268  }
269  }
270 
271  bool ret = num_wrong_dim_simplices == 0 && num_wrong_number_of_cofaces == 0;
272 
273  if (verbose_level >= 1) {
274  std::cerr << "Pure pseudo-manifold: ";
275  if (ret) {
276  std::cerr << green << "YES" << white << "\n";
277  } else {
278  std::cerr << red << "NO" << white << "\n"
279  << " * Number of wrong dimension simplices: "
280  << num_wrong_dim_simplices << "\n"
281  << " * Number of wrong number of cofaces: "
282  << num_wrong_number_of_cofaces << "\n";
283  }
284  }
285 
286  if (p_num_wrong_dim_simplices)
287  *p_num_wrong_dim_simplices = num_wrong_dim_simplices;
288  if (p_num_wrong_number_of_cofaces)
289  *p_num_wrong_number_of_cofaces = num_wrong_number_of_cofaces;
290 
291  return ret;
292  }
293 
294  template <int K>
295  std::size_t num_K_simplices() const {
296  Simplex_set k_simplices;
297 
298  for (Complex::const_iterator it_simplex = m_complex.begin(),
299  it_simplex_end = m_complex.end();
300  it_simplex != it_simplex_end;
301  ++it_simplex) {
302  if (it_simplex->size() == K + 1) {
303  k_simplices.insert(*it_simplex);
304  } else if (it_simplex->size() > K + 1) {
305  // Get the k-faces composing the simplex
306  combinations(
307  *it_simplex, K + 1, std::inserter(k_simplices, k_simplices.begin()));
308  }
309  }
310 
311  return k_simplices.size();
312  }
313 
314  std::ptrdiff_t euler_characteristic(bool verbose = false) const {
315  if (verbose)
316  std::cerr << "\nComputing Euler characteristic of the complex...\n";
317 
318  std::size_t num_vertices = num_K_simplices<0>();
319  std::size_t num_edges = num_K_simplices<1>();
320  std::size_t num_triangles = num_K_simplices<2>();
321 
322  std::ptrdiff_t ec =
323  (std::ptrdiff_t) num_vertices
324  - (std::ptrdiff_t) num_edges
325  + (std::ptrdiff_t) num_triangles;
326 
327  if (verbose)
328  std::cerr << "Euler characteristic: V - E + F = "
329  << num_vertices << " - " << num_edges << " + " << num_triangles << " = "
330  << blue
331  << ec
332  << white << "\n";
333 
334  return ec;
335  }
336 
337  // TODO(CJ): ADD COMMENTS
338 
339  bool is_pure_pseudomanifold(
340  int simplex_dim,
341  std::size_t num_vertices,
342  bool allow_borders = false,
343  bool exit_at_the_first_problem = false,
344  int verbose_level = 0,
345  std::size_t *p_num_wrong_dim_simplices = NULL,
346  std::size_t *p_num_wrong_number_of_cofaces = NULL,
347  std::size_t *p_num_unconnected_stars = NULL,
348  Simplex_set *p_wrong_dim_simplices = NULL,
349  Simplex_set *p_wrong_number_of_cofaces_simplices = NULL,
350  Simplex_set *p_unconnected_stars_simplices = NULL) const {
351  // If simplex_dim == 1, we do not need to check if stars are connected
352  if (simplex_dim == 1) {
353  if (p_num_unconnected_stars)
354  *p_num_unconnected_stars = 0;
355  return is_pure_pseudomanifold__do_not_check_if_stars_are_connected(simplex_dim,
356  allow_borders,
357  exit_at_the_first_problem,
358  verbose_level,
359  p_num_wrong_dim_simplices,
360  p_num_wrong_number_of_cofaces);
361  }
362  // Associates each vertex (= the index in the vector)
363  // to its star (list of simplices)
364  typedef std::vector<std::vector<Complex::const_iterator> > Stars;
365  std::size_t num_wrong_dim_simplices = 0;
366  std::size_t num_wrong_number_of_cofaces = 0;
367  std::size_t num_unconnected_stars = 0;
368 
369  // Fills a Stars data structure
370  Stars stars;
371  stars.resize(num_vertices);
372  for (Complex::const_iterator it_simplex = m_complex.begin(),
373  it_simplex_end = m_complex.end();
374  it_simplex != it_simplex_end;
375  ++it_simplex) {
376  if (static_cast<int> (it_simplex->size()) != simplex_dim + 1) {
377  if (verbose_level >= 2)
378  std::cerr << "Found a simplex with dim = "
379  << it_simplex->size() - 1 << "\n";
380  ++num_wrong_dim_simplices;
381  if (p_wrong_dim_simplices)
382  p_wrong_dim_simplices->insert(*it_simplex);
383  } else {
384  for (Simplex::const_iterator it_point_idx = it_simplex->begin();
385  it_point_idx != it_simplex->end();
386  ++it_point_idx) {
387  stars[*it_point_idx].push_back(it_simplex);
388  }
389  }
390  }
391 
392  // Now, for each star, we have a vector of its d-simplices
393  // i.e. one index for each d-simplex
394  // Boost Graph only deals with indexes, so we also need indexes for the
395  // (d-1)-simplices
396  std::size_t center_vertex_index = 0;
397  for (Stars::const_iterator it_star = stars.begin();
398  it_star != stars.end();
399  ++it_star, ++center_vertex_index) {
400  typedef std::map<Simplex, std::vector<std::size_t> >
401  Dm1_faces_to_adj_D_faces;
402  Dm1_faces_to_adj_D_faces dm1_faces_to_adj_d_faces;
403 
404  for (std::size_t i_dsimpl = 0; i_dsimpl < it_star->size(); ++i_dsimpl) {
405  Simplex dm1_simpl_of_link = *((*it_star)[i_dsimpl]);
406  dm1_simpl_of_link.erase(center_vertex_index);
407  // Copy it to a vector so that we can use operator[] on it
408  std::vector<std::size_t> dm1_simpl_of_link_vec(
409  dm1_simpl_of_link.begin(), dm1_simpl_of_link.end());
410 
411  CGAL::Combination_enumerator<int> dm2_simplices(
412  simplex_dim - 1, 0, simplex_dim);
413  for (; !dm2_simplices.finished(); ++dm2_simplices) {
414  Simplex dm2_simpl;
415  for (int j = 0; j < simplex_dim - 1; ++j)
416  dm2_simpl.insert(dm1_simpl_of_link_vec[dm2_simplices[j]]);
417  dm1_faces_to_adj_d_faces[dm2_simpl].push_back(i_dsimpl);
418  }
419  }
420 
421  Adj_graph adj_graph;
422  std::vector<Graph_vertex> d_faces_descriptors;
423  d_faces_descriptors.resize(it_star->size());
424  for (std::size_t j = 0; j < it_star->size(); ++j)
425  d_faces_descriptors[j] = boost::add_vertex(adj_graph);
426 
427  Dm1_faces_to_adj_D_faces::const_iterator dm1_to_d_it =
428  dm1_faces_to_adj_d_faces.begin();
429  Dm1_faces_to_adj_D_faces::const_iterator dm1_to_d_it_end =
430  dm1_faces_to_adj_d_faces.end();
431  for (std::size_t i_km1_face = 0;
432  dm1_to_d_it != dm1_to_d_it_end;
433  ++dm1_to_d_it, ++i_km1_face) {
434  Graph_vertex km1_gv = boost::add_vertex(adj_graph);
435 
436  for (std::vector<std::size_t>::const_iterator kface_it =
437  dm1_to_d_it->second.begin();
438  kface_it != dm1_to_d_it->second.end();
439  ++kface_it) {
440  boost::add_edge(km1_gv, *kface_it, adj_graph);
441  }
442 
443  if (dm1_to_d_it->second.size() != 2
444  && (!allow_borders || dm1_to_d_it->second.size() != 1)) {
445  ++num_wrong_number_of_cofaces;
446  if (p_wrong_number_of_cofaces_simplices) {
447  for (auto idx : dm1_to_d_it->second)
448  p_wrong_number_of_cofaces_simplices->insert(*((*it_star)[idx]));
449  }
450  }
451  }
452 
453  // What is left is to check the connexity
454  bool is_connected = true;
455  if (boost::num_vertices(adj_graph) > 0) {
456  std::vector<int> components(boost::num_vertices(adj_graph));
457  is_connected =
458  (boost::connected_components(adj_graph, &components[0]) == 1);
459  }
460 
461  if (!is_connected) {
462  if (verbose_level >= 2)
463  std::cerr << "Error: star #" << center_vertex_index
464  << " is not connected\n";
465  ++num_unconnected_stars;
466  if (p_unconnected_stars_simplices) {
467  for (std::vector<Complex::const_iterator>::const_iterator
468  it_simpl = it_star->begin(),
469  it_simpl_end = it_star->end();
470  it_simpl != it_simpl_end;
471  ++it_simpl) {
472  p_unconnected_stars_simplices->insert(**it_simpl);
473  }
474  }
475  }
476  }
477 
478  // Each one has been counted several times ("simplex_dim" times)
479  num_wrong_number_of_cofaces /= simplex_dim;
480 
481  bool ret =
482  num_wrong_dim_simplices == 0
483  && num_wrong_number_of_cofaces == 0
484  && num_unconnected_stars == 0;
485 
486  if (verbose_level >= 1) {
487  std::cerr << "Pure pseudo-manifold: ";
488  if (ret) {
489  std::cerr << green << "YES" << white << "\n";
490  } else {
491  std::cerr << red << "NO" << white << "\n"
492  << " * Number of wrong dimension simplices: "
493  << num_wrong_dim_simplices << "\n"
494  << " * Number of wrong number of cofaces: "
495  << num_wrong_number_of_cofaces << "\n"
496  << " * Number of not-connected stars: "
497  << num_unconnected_stars << "\n";
498  }
499  }
500 
501  if (p_num_wrong_dim_simplices)
502  *p_num_wrong_dim_simplices = num_wrong_dim_simplices;
503  if (p_num_wrong_number_of_cofaces)
504  *p_num_wrong_number_of_cofaces = num_wrong_number_of_cofaces;
505  if (p_num_unconnected_stars)
506  *p_num_unconnected_stars = num_unconnected_stars;
507 
508  return ret;
509  }
510 
511  private:
512  typedef Simplex_set Complex;
513 
514  // graph is an adjacency list
515  typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Adj_graph;
516  // map that gives to a certain simplex its node in graph and its dimension
517  typedef boost::graph_traits<Adj_graph>::vertex_descriptor Graph_vertex;
518  typedef boost::graph_traits<Adj_graph>::edge_descriptor Graph_edge;
519 
520  Complex m_complex;
521 }; // class Simplicial_complex
522 
523 } // namespace internal
524 } // namespace tangential_complex
525 } // namespace Gudhi
526 
527 #endif // TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_
GUDHI  Version 3.4.1  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Thu Dec 23 2021 15:37:05 for GUDHI by Doxygen 1.9.1