Simplicial_complex.h
1 /* This file is part of the Gudhi Library. The Gudhi library
2  * (Geometric Understanding in Higher Dimensions) is a generic C++
3  * library for computational topology.
4  *
5  * Author(s): Clement Jamin
6  *
7  * Copyright (C) 2016 Inria
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with this program. If not, see <http://www.gnu.org/licenses/>.
21  */
22 
23 #ifndef TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_
24 #define TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_
25 
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>
30 
31 #include <CGAL/iterator.h>
32 
33 // For is_pure_pseudomanifold
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>
38 
39 #include <algorithm>
40 #include <string>
41 #include <fstream>
42 #include <map> // for map<>
43 #include <vector> // for vector<>
44 #include <set> // for set<>
45 
46 namespace Gudhi {
47 namespace tangential_complex {
48 namespace internal {
49 
50 class Simplicial_complex {
51  public:
52  typedef boost::container::flat_set<std::size_t> Simplex;
53  typedef std::set<Simplex> Simplex_set;
54 
55  // If perform_checks = true, the function:
56  // - won't insert the simplex if it is already in a higher dim simplex
57  // - will erase any lower-dim simplices that are faces of the new simplex
58  // Returns true if the simplex was added
59  bool add_simplex(
60  const Simplex &s, bool perform_checks = true) {
61  if (perform_checks) {
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;
68  ++it_simplex) {
69  // Check if the simplex is not already in a higher dim simplex
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())) {
74  // No need to insert it, then
75  return false;
76  }
77  // Check if the simplex includes some lower-dim simplices
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);
82  // We don't need to check higher-sim simplices any more
83  check_higher_dim_simpl = false;
84  }
85  }
86  for (std::vector<Complex::iterator>::const_iterator it = to_erase.begin();
87  it != to_erase.end(); ++it) {
88  m_complex.erase(*it);
89  }
90  }
91  return m_complex.insert(s).second;
92  }
93 
94  const Simplex_set &simplex_range() const {
95  return m_complex;
96  }
97 
98  bool empty() {
99  return m_complex.empty();
100  }
101 
102  void clear() {
103  m_complex.clear();
104  }
105 
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;
111  ++it_simplex) {
112  if (test(*it_simplex))
113  *out++ = *it_simplex;
114  }
115  }
116 
117  // When a simplex S has only one co-face C, we can remove S and C
118  // without changing the topology
119 
120  void collapse(int max_simplex_dim, bool quiet = false) {
121 #ifdef DEBUG_TRACES
122  if (!quiet)
123  std::cerr << "Collapsing... ";
124 #endif
125  // We note k = max_simplex_dim - 1
126  int k = max_simplex_dim - 1;
127 
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;
131 
132  std::size_t num_collapsed_maximal_simplices = 0;
133  do {
134  num_collapsed_maximal_simplices = 0;
135  // Create a map associating each non-maximal k-faces to the list of its
136  // maximal cofaces
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;
141  ++it_simplex) {
142  if (static_cast<int> (it_simplex->size()) > k + 1) {
143  std::vector<Simplex> k_faces;
144  // Get the k-faces composing the simplex
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);
148  }
149  }
150 
151  // For each non-maximal k-face F, if F has only one maximal coface Cf:
152  // - Look for the other k-faces F2, F3... of Cf in the map and:
153  // * if the list contains only Cf, clear the list (we don't remove the
154  // list since it creates troubles with the iterators) and add the F2,
155  // F3... to the complex
156  // * otherwise, remove Cf from the associated list
157  // - Remove Cf from the complex
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;
161  ++it_map_elt) {
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"));
168  // Get the k-faces composing the simplex
169  combinations(*it_Cf, k + 1, std::back_inserter(k_faces));
170  for (const auto &f2 : k_faces) {
171  // Skip F
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);
177  } else { // it_comb_in_map->second.size() > 1
178  Simplex_iterator_list::iterator it = std::find(it_comb_in_map->second.begin(),
179  it_comb_in_map->second.end(),
180  it_Cf);
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);
184  }
185  }
186  }
187  m_complex.erase(it_Cf);
188  ++num_collapsed_maximal_simplices;
189  }
190  }
191  // Repeat until no maximal simplex got removed
192  } while (num_collapsed_maximal_simplices > 0);
193 
194  // Collapse the lower dimension simplices
195  if (k > 0)
196  collapse(max_simplex_dim - 1, true);
197 
198 #ifdef DEBUG_TRACES
199  if (!quiet)
200  std::cerr << "done.\n";
201 #endif
202  }
203 
204  void display_stats() const {
205  std::cerr << yellow << "Complex stats:\n" << white;
206 
207  if (m_complex.empty()) {
208  std::cerr << " * No simplices.\n";
209  } else {
210  // Number of simplex for each dimension
211  std::map<int, std::size_t> simplex_stats;
212 
213  for (Complex::const_iterator it_simplex = m_complex.begin(),
214  it_simplex_end = m_complex.end();
215  it_simplex != it_simplex_end;
216  ++it_simplex) {
217  ++simplex_stats[static_cast<int> (it_simplex->size()) - 1];
218  }
219 
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";
224  }
225  }
226  }
227 
228  // verbose_level = 0, 1 or 2
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;
237 
238  std::size_t num_wrong_dim_simplices = 0;
239  std::size_t num_wrong_number_of_cofaces = 0;
240 
241  // Counts the number of cofaces of each K_1_face
242 
243  // Create a map associating each non-maximal k-faces to the list of its
244  // maximal cofaces
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;
249  ++it_simplex) {
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;
255  } else {
256  std::vector<K_1_face> k_1_faces;
257  // Get the facets composing the simplex
258  combinations(
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];
262  }
263  }
264  }
265 
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;
269  ++it_map_elt) {
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";
275 
276  if (exit_at_the_first_problem)
277  return false;
278  else
279  ++num_wrong_number_of_cofaces;
280  }
281  }
282 
283  bool ret = num_wrong_dim_simplices == 0 && num_wrong_number_of_cofaces == 0;
284 
285  if (verbose_level >= 1) {
286  std::cerr << "Pure pseudo-manifold: ";
287  if (ret) {
288  std::cerr << green << "YES" << white << "\n";
289  } else {
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";
295  }
296  }
297 
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;
302 
303  return ret;
304  }
305 
306  template <int K>
307  std::size_t num_K_simplices() const {
308  Simplex_set k_simplices;
309 
310  for (Complex::const_iterator it_simplex = m_complex.begin(),
311  it_simplex_end = m_complex.end();
312  it_simplex != it_simplex_end;
313  ++it_simplex) {
314  if (it_simplex->size() == K + 1) {
315  k_simplices.insert(*it_simplex);
316  } else if (it_simplex->size() > K + 1) {
317  // Get the k-faces composing the simplex
318  combinations(
319  *it_simplex, K + 1, std::inserter(k_simplices, k_simplices.begin()));
320  }
321  }
322 
323  return k_simplices.size();
324  }
325 
326  std::ptrdiff_t euler_characteristic(bool verbose = false) const {
327  if (verbose)
328  std::cerr << "\nComputing Euler characteristic of the complex...\n";
329 
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>();
333 
334  std::ptrdiff_t ec =
335  (std::ptrdiff_t) num_vertices
336  - (std::ptrdiff_t) num_edges
337  + (std::ptrdiff_t) num_triangles;
338 
339  if (verbose)
340  std::cerr << "Euler characteristic: V - E + F = "
341  << num_vertices << " - " << num_edges << " + " << num_triangles << " = "
342  << blue
343  << ec
344  << white << "\n";
345 
346  return ec;
347  }
348 
349  // TODO(CJ): ADD COMMENTS
350 
351  bool is_pure_pseudomanifold(
352  int simplex_dim,
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 {
363  // If simplex_dim == 1, we do not need to check if stars are connected
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,
368  allow_borders,
369  exit_at_the_first_problem,
370  verbose_level,
371  p_num_wrong_dim_simplices,
372  p_num_wrong_number_of_cofaces);
373  }
374  // Associates each vertex (= the index in the vector)
375  // to its star (list of simplices)
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;
380 
381  // Fills a Stars data structure
382  Stars stars;
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;
387  ++it_simplex) {
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);
395  } else {
396  for (Simplex::const_iterator it_point_idx = it_simplex->begin();
397  it_point_idx != it_simplex->end();
398  ++it_point_idx) {
399  stars[*it_point_idx].push_back(it_simplex);
400  }
401  }
402  }
403 
404  // Now, for each star, we have a vector of its d-simplices
405  // i.e. one index for each d-simplex
406  // Boost Graph only deals with indexes, so we also need indexes for the
407  // (d-1)-simplices
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;
415 
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);
419  // Copy it to a vector so that we can use operator[] on it
420  std::vector<std::size_t> dm1_simpl_of_link_vec(
421  dm1_simpl_of_link.begin(), dm1_simpl_of_link.end());
422 
423  CGAL::Combination_enumerator<int> dm2_simplices(
424  simplex_dim - 1, 0, simplex_dim);
425  for (; !dm2_simplices.finished(); ++dm2_simplices) {
426  Simplex dm2_simpl;
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);
430  }
431  }
432 
433  Adj_graph adj_graph;
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);
438 
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);
447 
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();
451  ++kface_it) {
452  boost::add_edge(km1_gv, *kface_it, adj_graph);
453  }
454 
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]));
461  }
462  }
463  }
464 
465  // What is left is to check the connexity
466  bool is_connected = true;
467  if (boost::num_vertices(adj_graph) > 0) {
468  std::vector<int> components(boost::num_vertices(adj_graph));
469  is_connected =
470  (boost::connected_components(adj_graph, &components[0]) == 1);
471  }
472 
473  if (!is_connected) {
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;
483  ++it_simpl) {
484  p_unconnected_stars_simplices->insert(**it_simpl);
485  }
486  }
487  }
488  }
489 
490  // Each one has been counted several times ("simplex_dim" times)
491  num_wrong_number_of_cofaces /= simplex_dim;
492 
493  bool ret =
494  num_wrong_dim_simplices == 0
495  && num_wrong_number_of_cofaces == 0
496  && num_unconnected_stars == 0;
497 
498  if (verbose_level >= 1) {
499  std::cerr << "Pure pseudo-manifold: ";
500  if (ret) {
501  std::cerr << green << "YES" << white << "\n";
502  } else {
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";
510  }
511  }
512 
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;
519 
520  return ret;
521  }
522 
523  private:
524  typedef Simplex_set Complex;
525 
526  // graph is an adjacency list
527  typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Adj_graph;
528  // map that gives to a certain simplex its node in graph and its dimension
529  typedef boost::graph_traits<Adj_graph>::vertex_descriptor Graph_vertex;
530  typedef boost::graph_traits<Adj_graph>::edge_descriptor Graph_edge;
531 
532  Complex m_complex;
533 }; // class Simplicial_complex
534 
535 } // namespace internal
536 } // namespace tangential_complex
537 } // namespace Gudhi
538 
539 #endif // TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_
Definition: SimplicialComplexForAlpha.h:26
GUDHI  Version 2.3.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : GPL v3 Generated on Fri Jan 22 2021 03:21:44 for GUDHI by Doxygen 1.8.13