Loading...
Searching...
No Matches
Simplicial_complex.h
/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
* See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
* Author(s): Clement Jamin
*
* Copyright (C) 2016 Inria
*
* Modification(s):
* - YYYY/MM Author: Description of the modification
*/
#ifndef TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_
#define TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_
#include <gudhi/Tangential_complex/config.h>
#include <gudhi/Tangential_complex/utilities.h>
#include <gudhi/Debug_utils.h>
#include <gudhi/console_color.h>
#include <CGAL/iterator.h>
// For is_pure_pseudomanifold
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/connected_components.hpp>
#include <boost/container/flat_set.hpp>
#include <algorithm>
#include <string>
#include <fstream>
#include <map> // for map<>
#include <vector> // for vector<>
#include <set> // for set<>
namespace Gudhi {
namespace tangential_complex {
namespace internal {
class Simplicial_complex {
public:
typedef boost::container::flat_set<std::size_t> Simplex;
typedef std::set<Simplex> Simplex_set;
// If perform_checks = true, the function:
// - won't insert the simplex if it is already in a higher dim simplex
// - will erase any lower-dim simplices that are faces of the new simplex
// Returns true if the simplex was added
bool add_simplex(
const Simplex &s, bool perform_checks = true) {
if (perform_checks) {
unsigned int num_pts = static_cast<int> (s.size());
std::vector<Complex::iterator> to_erase;
bool check_higher_dim_simpl = true;
for (Complex::iterator it_simplex = m_complex.begin(),
it_simplex_end = m_complex.end();
it_simplex != it_simplex_end;
++it_simplex) {
// Check if the simplex is not already in a higher dim simplex
if (check_higher_dim_simpl
&& it_simplex->size() > num_pts
&& std::includes(it_simplex->begin(), it_simplex->end(),
s.begin(), s.end())) {
// No need to insert it, then
return false;
}
// Check if the simplex includes some lower-dim simplices
if (it_simplex->size() < num_pts
&& std::includes(s.begin(), s.end(),
it_simplex->begin(), it_simplex->end())) {
to_erase.push_back(it_simplex);
// We don't need to check higher-sim simplices any more
check_higher_dim_simpl = false;
}
}
for (std::vector<Complex::iterator>::const_iterator it = to_erase.begin();
it != to_erase.end(); ++it) {
m_complex.erase(*it);
}
}
return m_complex.insert(s).second;
}
const Simplex_set &simplex_range() const {
return m_complex;
}
bool empty() {
return m_complex.empty();
}
void clear() {
m_complex.clear();
}
template <typename Test, typename Output_it>
void get_simplices_matching_test(Test test, Output_it out) {
for (Complex::const_iterator it_simplex = m_complex.begin(),
it_simplex_end = m_complex.end();
it_simplex != it_simplex_end;
++it_simplex) {
if (test(*it_simplex))
*out++ = *it_simplex;
}
}
// When a simplex S has only one co-face C, we can remove S and C
// without changing the topology
void collapse(int max_simplex_dim, bool quiet = false) {
#ifdef DEBUG_TRACES
if (!quiet)
std::cerr << "Collapsing... ";
#endif
// We note k = max_simplex_dim - 1
int k = max_simplex_dim - 1;
typedef Complex::iterator Simplex_iterator;
typedef std::vector<Simplex_iterator> Simplex_iterator_list;
typedef std::map<Simplex, Simplex_iterator_list> Cofaces_map;
std::size_t num_collapsed_maximal_simplices = 0;
do {
num_collapsed_maximal_simplices = 0;
// Create a map associating each non-maximal k-faces to the list of its
// maximal cofaces
Cofaces_map cofaces_map;
for (Complex::const_iterator it_simplex = m_complex.begin(),
it_simplex_end = m_complex.end();
it_simplex != it_simplex_end;
++it_simplex) {
if (static_cast<int> (it_simplex->size()) > k + 1) {
std::vector<Simplex> k_faces;
// Get the k-faces composing the simplex
combinations(*it_simplex, k + 1, std::back_inserter(k_faces));
for (const auto &comb : k_faces)
cofaces_map[comb].push_back(it_simplex);
}
}
// For each non-maximal k-face F, if F has only one maximal coface Cf:
// - Look for the other k-faces F2, F3... of Cf in the map and:
// * if the list contains only Cf, clear the list (we don't remove the
// list since it creates troubles with the iterators) and add the F2,
// F3... to the complex
// * otherwise, remove Cf from the associated list
// - Remove Cf from the complex
for (Cofaces_map::const_iterator it_map_elt = cofaces_map.begin(),
it_map_end = cofaces_map.end();
it_map_elt != it_map_end;
++it_map_elt) {
if (it_map_elt->second.size() == 1) {
std::vector<Simplex> k_faces;
const Simplex_iterator_list::value_type &it_Cf =
*it_map_elt->second.begin();
GUDHI_CHECK(it_Cf->size() == max_simplex_dim + 1,
std::logic_error("Wrong dimension"));
// Get the k-faces composing the simplex
combinations(*it_Cf, k + 1, std::back_inserter(k_faces));
for (const auto &f2 : k_faces) {
// Skip F
if (f2 != it_map_elt->first) {
Cofaces_map::iterator it_comb_in_map = cofaces_map.find(f2);
if (it_comb_in_map->second.size() == 1) {
it_comb_in_map->second.clear();
m_complex.insert(f2);
} else { // it_comb_in_map->second.size() > 1
Simplex_iterator_list::iterator it = std::find(it_comb_in_map->second.begin(),
it_comb_in_map->second.end(),
it_Cf);
GUDHI_CHECK(it != it_comb_in_map->second.end(),
std::logic_error("Error: it == it_comb_in_map->second.end()"));
it_comb_in_map->second.erase(it);
}
}
}
m_complex.erase(it_Cf);
++num_collapsed_maximal_simplices;
}
}
// Repeat until no maximal simplex got removed
} while (num_collapsed_maximal_simplices > 0);
// Collapse the lower dimension simplices
if (k > 0)
collapse(max_simplex_dim - 1, true);
#ifdef DEBUG_TRACES
if (!quiet)
std::cerr << "done.\n";
#endif
}
void display_stats() const {
std::cerr << yellow << "Complex stats:\n" << white;
if (m_complex.empty()) {
std::cerr << " * No simplices.\n";
} else {
// Number of simplex for each dimension
std::map<int, std::size_t> simplex_stats;
for (Complex::const_iterator it_simplex = m_complex.begin(),
it_simplex_end = m_complex.end();
it_simplex != it_simplex_end;
++it_simplex) {
++simplex_stats[static_cast<int> (it_simplex->size()) - 1];
}
for (std::map<int, std::size_t>::const_iterator it_map = simplex_stats.begin();
it_map != simplex_stats.end(); ++it_map) {
std::cerr << " * " << it_map->first << "-simplices: "
<< it_map->second << "\n";
}
}
}
// verbose_level = 0, 1 or 2
bool is_pure_pseudomanifold__do_not_check_if_stars_are_connected(int simplex_dim,
bool allow_borders = false,
bool exit_at_the_first_problem = false,
int verbose_level = 0,
std::size_t *p_num_wrong_dim_simplices = NULL,
std::size_t *p_num_wrong_number_of_cofaces = NULL) const {
typedef Simplex K_1_face;
typedef std::map<K_1_face, std::size_t> Cofaces_map;
std::size_t num_wrong_dim_simplices = 0;
std::size_t num_wrong_number_of_cofaces = 0;
// Counts the number of cofaces of each K_1_face
// Create a map associating each non-maximal k-faces to the list of its
// maximal cofaces
Cofaces_map cofaces_map;
for (Complex::const_iterator it_simplex = m_complex.begin(),
it_simplex_end = m_complex.end();
it_simplex != it_simplex_end;
++it_simplex) {
if (static_cast<int> (it_simplex->size()) != simplex_dim + 1) {
if (verbose_level >= 2)
std::cerr << "Found a simplex with dim = "
<< it_simplex->size() - 1 << "\n";
++num_wrong_dim_simplices;
} else {
std::vector<K_1_face> k_1_faces;
// Get the facets composing the simplex
combinations(
*it_simplex, simplex_dim, std::back_inserter(k_1_faces));
for (const auto &k_1_face : k_1_faces) {
++cofaces_map[k_1_face];
}
}
}
for (Cofaces_map::const_iterator it_map_elt = cofaces_map.begin(),
it_map_end = cofaces_map.end();
it_map_elt != it_map_end;
++it_map_elt) {
if (it_map_elt->second != 2
&& (!allow_borders || it_map_elt->second != 1)) {
if (verbose_level >= 2)
std::cerr << "Found a k-1-face with "
<< it_map_elt->second << " cofaces\n";
if (exit_at_the_first_problem)
return false;
else
++num_wrong_number_of_cofaces;
}
}
bool ret = num_wrong_dim_simplices == 0 && num_wrong_number_of_cofaces == 0;
if (verbose_level >= 1) {
std::cerr << "Pure pseudo-manifold: ";
if (ret) {
std::cerr << green << "YES" << white << "\n";
} else {
std::cerr << red << "NO" << white << "\n"
<< " * Number of wrong dimension simplices: "
<< num_wrong_dim_simplices << "\n"
<< " * Number of wrong number of cofaces: "
<< num_wrong_number_of_cofaces << "\n";
}
}
if (p_num_wrong_dim_simplices)
*p_num_wrong_dim_simplices = num_wrong_dim_simplices;
if (p_num_wrong_number_of_cofaces)
*p_num_wrong_number_of_cofaces = num_wrong_number_of_cofaces;
return ret;
}
template <int K>
std::size_t num_K_simplices() const {
Simplex_set k_simplices;
for (Complex::const_iterator it_simplex = m_complex.begin(),
it_simplex_end = m_complex.end();
it_simplex != it_simplex_end;
++it_simplex) {
if (it_simplex->size() == K + 1) {
k_simplices.insert(*it_simplex);
} else if (it_simplex->size() > K + 1) {
// Get the k-faces composing the simplex
combinations(
*it_simplex, K + 1, std::inserter(k_simplices, k_simplices.begin()));
}
}
return k_simplices.size();
}
std::ptrdiff_t euler_characteristic(bool verbose = false) const {
if (verbose)
std::cerr << "\nComputing Euler characteristic of the complex...\n";
std::size_t num_vertices = num_K_simplices<0>();
std::size_t num_edges = num_K_simplices<1>();
std::size_t num_triangles = num_K_simplices<2>();
std::ptrdiff_t ec =
(std::ptrdiff_t) num_vertices
- (std::ptrdiff_t) num_edges
+ (std::ptrdiff_t) num_triangles;
if (verbose)
std::cerr << "Euler characteristic: V - E + F = "
<< num_vertices << " - " << num_edges << " + " << num_triangles << " = "
<< blue
<< ec
<< white << "\n";
return ec;
}
// TODO(CJ): ADD COMMENTS
bool is_pure_pseudomanifold(
int simplex_dim,
std::size_t num_vertices,
bool allow_borders = false,
bool exit_at_the_first_problem = false,
int verbose_level = 0,
std::size_t *p_num_wrong_dim_simplices = NULL,
std::size_t *p_num_wrong_number_of_cofaces = NULL,
std::size_t *p_num_unconnected_stars = NULL,
Simplex_set *p_wrong_dim_simplices = NULL,
Simplex_set *p_wrong_number_of_cofaces_simplices = NULL,
Simplex_set *p_unconnected_stars_simplices = NULL) const {
// If simplex_dim == 1, we do not need to check if stars are connected
if (simplex_dim == 1) {
if (p_num_unconnected_stars)
*p_num_unconnected_stars = 0;
return is_pure_pseudomanifold__do_not_check_if_stars_are_connected(simplex_dim,
allow_borders,
exit_at_the_first_problem,
verbose_level,
p_num_wrong_dim_simplices,
p_num_wrong_number_of_cofaces);
}
// Associates each vertex (= the index in the vector)
// to its star (list of simplices)
typedef std::vector<std::vector<Complex::const_iterator> > Stars;
std::size_t num_wrong_dim_simplices = 0;
std::size_t num_wrong_number_of_cofaces = 0;
std::size_t num_unconnected_stars = 0;
// Fills a Stars data structure
Stars stars;
stars.resize(num_vertices);
for (Complex::const_iterator it_simplex = m_complex.begin(),
it_simplex_end = m_complex.end();
it_simplex != it_simplex_end;
++it_simplex) {
if (static_cast<int> (it_simplex->size()) != simplex_dim + 1) {
if (verbose_level >= 2)
std::cerr << "Found a simplex with dim = "
<< it_simplex->size() - 1 << "\n";
++num_wrong_dim_simplices;
if (p_wrong_dim_simplices)
p_wrong_dim_simplices->insert(*it_simplex);
} else {
for (Simplex::const_iterator it_point_idx = it_simplex->begin();
it_point_idx != it_simplex->end();
++it_point_idx) {
stars[*it_point_idx].push_back(it_simplex);
}
}
}
// Now, for each star, we have a vector of its d-simplices
// i.e. one index for each d-simplex
// Boost Graph only deals with indexes, so we also need indexes for the
// (d-1)-simplices
std::size_t center_vertex_index = 0;
for (Stars::const_iterator it_star = stars.begin();
it_star != stars.end();
++it_star, ++center_vertex_index) {
typedef std::map<Simplex, std::vector<std::size_t> >
Dm1_faces_to_adj_D_faces;
Dm1_faces_to_adj_D_faces dm1_faces_to_adj_d_faces;
for (std::size_t i_dsimpl = 0; i_dsimpl < it_star->size(); ++i_dsimpl) {
Simplex dm1_simpl_of_link = *((*it_star)[i_dsimpl]);
dm1_simpl_of_link.erase(center_vertex_index);
// Copy it to a vector so that we can use operator[] on it
std::vector<std::size_t> dm1_simpl_of_link_vec(
dm1_simpl_of_link.begin(), dm1_simpl_of_link.end());
CGAL::Combination_enumerator<int> dm2_simplices(
simplex_dim - 1, 0, simplex_dim);
for (; !dm2_simplices.finished(); ++dm2_simplices) {
Simplex dm2_simpl;
for (int j = 0; j < simplex_dim - 1; ++j)
dm2_simpl.insert(dm1_simpl_of_link_vec[dm2_simplices[j]]);
dm1_faces_to_adj_d_faces[dm2_simpl].push_back(i_dsimpl);
}
}
Adj_graph adj_graph;
std::vector<Graph_vertex> d_faces_descriptors;
d_faces_descriptors.resize(it_star->size());
for (std::size_t j = 0; j < it_star->size(); ++j)
d_faces_descriptors[j] = boost::add_vertex(adj_graph);
Dm1_faces_to_adj_D_faces::const_iterator dm1_to_d_it =
dm1_faces_to_adj_d_faces.begin();
Dm1_faces_to_adj_D_faces::const_iterator dm1_to_d_it_end =
dm1_faces_to_adj_d_faces.end();
for (std::size_t i_km1_face = 0;
dm1_to_d_it != dm1_to_d_it_end;
++dm1_to_d_it, ++i_km1_face) {
Graph_vertex km1_gv = boost::add_vertex(adj_graph);
for (std::vector<std::size_t>::const_iterator kface_it =
dm1_to_d_it->second.begin();
kface_it != dm1_to_d_it->second.end();
++kface_it) {
boost::add_edge(km1_gv, *kface_it, adj_graph);
}
if (dm1_to_d_it->second.size() != 2
&& (!allow_borders || dm1_to_d_it->second.size() != 1)) {
++num_wrong_number_of_cofaces;
if (p_wrong_number_of_cofaces_simplices) {
for (auto idx : dm1_to_d_it->second)
p_wrong_number_of_cofaces_simplices->insert(*((*it_star)[idx]));
}
}
}
// What is left is to check the connexity
bool is_connected = true;
if (boost::num_vertices(adj_graph) > 0) {
std::vector<int> components(boost::num_vertices(adj_graph));
is_connected =
(boost::connected_components(adj_graph, &components[0]) == 1);
}
if (!is_connected) {
if (verbose_level >= 2)
std::cerr << "Error: star #" << center_vertex_index
<< " is not connected\n";
++num_unconnected_stars;
if (p_unconnected_stars_simplices) {
for (std::vector<Complex::const_iterator>::const_iterator
it_simpl = it_star->begin(),
it_simpl_end = it_star->end();
it_simpl != it_simpl_end;
++it_simpl) {
p_unconnected_stars_simplices->insert(**it_simpl);
}
}
}
}
// Each one has been counted several times ("simplex_dim" times)
num_wrong_number_of_cofaces /= simplex_dim;
bool ret =
num_wrong_dim_simplices == 0
&& num_wrong_number_of_cofaces == 0
&& num_unconnected_stars == 0;
if (verbose_level >= 1) {
std::cerr << "Pure pseudo-manifold: ";
if (ret) {
std::cerr << green << "YES" << white << "\n";
} else {
std::cerr << red << "NO" << white << "\n"
<< " * Number of wrong dimension simplices: "
<< num_wrong_dim_simplices << "\n"
<< " * Number of wrong number of cofaces: "
<< num_wrong_number_of_cofaces << "\n"
<< " * Number of not-connected stars: "
<< num_unconnected_stars << "\n";
}
}
if (p_num_wrong_dim_simplices)
*p_num_wrong_dim_simplices = num_wrong_dim_simplices;
if (p_num_wrong_number_of_cofaces)
*p_num_wrong_number_of_cofaces = num_wrong_number_of_cofaces;
if (p_num_unconnected_stars)
*p_num_unconnected_stars = num_unconnected_stars;
return ret;
}
private:
typedef Simplex_set Complex;
// graph is an adjacency list
typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Adj_graph;
// map that gives to a certain simplex its node in graph and its dimension
typedef boost::graph_traits<Adj_graph>::vertex_descriptor Graph_vertex;
typedef boost::graph_traits<Adj_graph>::edge_descriptor Graph_edge;
Complex m_complex;
}; // class Simplicial_complex
} // namespace internal
} // namespace tangential_complex
} // namespace Gudhi
#endif // TANGENTIAL_COMPLEX_SIMPLICIAL_COMPLEX_H_
Gudhi
Gudhi namespace.
Definition SimplicialComplexForAlpha.h:14