11 #ifndef SIMPLEX_TREE_H_
12 #define SIMPLEX_TREE_H_
14 #include <gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h>
15 #include <gudhi/Simplex_tree/Simplex_tree_siblings.h>
16 #include <gudhi/Simplex_tree/Simplex_tree_iterators.h>
17 #include <gudhi/Simplex_tree/indexing_tag.h>
20 #include <gudhi/graph_simplicial_complex.h>
21 #include <gudhi/Debug_utils.h>
23 #include <boost/container/flat_map.hpp>
24 #include <boost/iterator/transform_iterator.hpp>
25 #include <boost/graph/adjacency_list.hpp>
26 #include <boost/range/adaptor/reversed.hpp>
29 #include <tbb/parallel_sort.h>
37 #include <initializer_list>
44 struct Simplex_tree_options_full_featured;
59 template<
typename SimplexTreeOptions = Simplex_tree_options_full_featured>
78 typedef Simplex_tree_node_explicit_storage<Simplex_tree> Node;
83 typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary;
87 typedef Simplex_tree_siblings<Simplex_tree, Dictionary> Siblings;
89 struct Key_simplex_base_real {
90 Key_simplex_base_real() : key_(-1) {}
96 struct Key_simplex_base_dummy {
97 Key_simplex_base_dummy() {}
102 typedef typename std::conditional<Options::store_key, Key_simplex_base_real, Key_simplex_base_dummy>::type
105 struct Filtration_simplex_base_real {
106 Filtration_simplex_base_real() : filt_(0) {}
112 struct Filtration_simplex_base_dummy {
113 Filtration_simplex_base_dummy() {}
118 Filtration_simplex_base_dummy>::type Filtration_simplex_base;
126 typedef typename Dictionary::iterator Dictionary_it;
127 typedef typename Dictionary_it::value_type Dit_value_t;
129 struct return_first {
193 boost::make_transform_iterator(root_.members_.begin(), return_first()),
194 boost::make_transform_iterator(root_.members_.end(), return_first()));
237 if (filtration_vect_.empty()) {
240 return filtration_vect_;
269 template<
class SimplexHandle>
282 root_(nullptr, null_vertex_),
289 std::cout <<
"Simplex_tree copy constructor" << std::endl;
290 #endif // DEBUG_TRACES
291 copy_from(complex_source);
299 std::cout <<
"Simplex_tree move constructor" << std::endl;
300 #endif // DEBUG_TRACES
301 move_from(complex_source);
305 complex_source.dimension_ = -1;
310 root_members_recursive_deletion();
316 std::cout <<
"Simplex_tree copy assignment" << std::endl;
317 #endif // DEBUG_TRACES
319 if (&complex_source !=
this) {
321 root_members_recursive_deletion();
323 copy_from(complex_source);
333 std::cout <<
"Simplex_tree move assignment" << std::endl;
334 #endif // DEBUG_TRACES
336 if (&complex_source !=
this) {
338 root_members_recursive_deletion();
340 move_from(complex_source);
349 null_vertex_ = complex_source.null_vertex_;
350 filtration_vect_.clear();
351 dimension_ = complex_source.dimension_;
352 auto root_source = complex_source.root_;
355 root_.members() = Dictionary(boost::container::ordered_unique_range, root_source.members().begin(), root_source.members().end());
357 for (
auto& map_el : root_.members()) {
358 map_el.second.assign_children(&root_);
360 rec_copy(&root_, &root_source);
364 void rec_copy(Siblings *sib, Siblings *sib_source) {
365 for (
auto sh = sib->members().begin(), sh_source = sib_source->members().begin();
366 sh != sib->members().end(); ++sh, ++sh_source) {
368 Siblings * newsib =
new Siblings(sib, sh_source->first);
369 newsib->members_.reserve(sh_source->second.children()->members().size());
370 for (
auto & child : sh_source->second.children()->members())
371 newsib->members_.emplace_hint(newsib->members_.end(), child.first, Node(newsib, child.second.filtration()));
372 rec_copy(newsib, sh_source->second.children());
373 sh->second.assign_children(newsib);
380 null_vertex_ = std::move(complex_source.null_vertex_);
381 root_ = std::move(complex_source.root_);
382 filtration_vect_ = std::move(complex_source.filtration_vect_);
383 dimension_ = std::move(complex_source.dimension_);
386 for (
auto& map_el : root_.members()) {
387 if (map_el.second.children() != &(complex_source.root_)) {
389 map_el.second.children()->oncles_ = &root_;
392 GUDHI_CHECK(map_el.second.children()->oncles_ ==
nullptr,
393 std::invalid_argument(
"Simplex_tree move constructor from an invalid Simplex_tree"));
395 map_el.second.assign_children(&root_);
401 void root_members_recursive_deletion() {
402 for (
auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) {
404 rec_delete(sh->second.children());
407 root_.members().clear();
411 void rec_delete(Siblings * sib) {
412 for (
auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
414 rec_delete(sh->second.children());
423 if ((null_vertex_ != st2.null_vertex_) ||
424 (dimension_ != st2.dimension_))
426 return rec_equal(&root_, &st2.root_);
431 return (!(*
this == st2));
436 bool rec_equal(Siblings* s1, Siblings* s2) {
437 if (s1->members().size() != s2->members().size())
439 for (
auto sh1 = s1->members().begin(), sh2 = s2->members().begin();
440 (sh1 != s1->members().end() && sh2 != s2->members().end()); ++sh1, ++sh2) {
441 if (sh1->first != sh2->first || sh1->second.filtration() != sh2->second.filtration())
447 if (!rec_equal(sh1->second.children(), sh2->second.children()))
460 return sh->second.key();
469 return filtration_vect_[idx];
479 return sh->second.filtration();
481 return std::numeric_limits<Filtration_value>::infinity();
490 std::invalid_argument(
"Simplex_tree::assign_filtration - cannot assign filtration on null_simplex"));
491 sh->second.assign_filtration(fv);
499 return Dictionary_it(
nullptr);
516 return root_.members_.size();
528 auto sib_begin = sib->members().begin();
529 auto sib_end = sib->members().end();
530 size_t simplices_number = sib_end - sib_begin;
531 for (
auto sh = sib_begin; sh != sib_end; ++sh) {
536 return simplices_number;
546 while (curr_sib !=
nullptr) {
548 curr_sib = curr_sib->oncles();
563 if (dimension_to_be_lowered_)
564 lower_upper_bound_dimension();
570 template<
class SimplexHandle>
573 return (sh->second.children()->parent() == sh->first);
583 template<
class InputVertexRange = std::initializer_list<Vertex_handle>>
585 auto first = std::begin(s);
586 auto last = std::end(s);
592 std::vector<Vertex_handle> copy(first, last);
593 std::sort(std::begin(copy), std::end(copy));
594 return find_simplex(copy);
600 Siblings * tmp_sib = &root_;
601 Dictionary_it tmp_dit;
605 GUDHI_CHECK(contiguous_vertices(),
"non-contiguous vertices");
607 if(v < 0 || v >= static_cast<Vertex_handle>(root_.members_.size()))
609 tmp_dit = root_.members_.begin() + v;
614 tmp_sib = tmp_dit->second.children();
617 tmp_dit = tmp_sib->members_.find(*vi++);
618 if (tmp_dit == tmp_sib->members_.end())
624 tmp_sib = tmp_dit->second.children();
632 assert(contiguous_vertices());
633 return root_.members_.begin() + v;
635 return root_.members_.find(v);
641 bool contiguous_vertices()
const {
642 if (root_.members_.empty())
return true;
643 if (root_.members_.begin()->first != 0)
return false;
644 if (std::prev(root_.members_.end())->first != static_cast<Vertex_handle>(root_.members_.size() - 1))
return false;
663 std::pair<Simplex_handle, bool> insert_vertex_vector(
const std::vector<Vertex_handle>&
simplex,
665 Siblings * curr_sib = &root_;
666 std::pair<Simplex_handle, bool> res_insert;
668 for (; vi !=
simplex.end() - 1; ++vi) {
669 GUDHI_CHECK(*vi !=
null_vertex(),
"cannot use the dummy null_vertex() as a real vertex");
670 res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib,
filtration));
672 res_insert.first->second.assign_children(
new Siblings(curr_sib, *vi));
674 curr_sib = res_insert.first->second.children();
676 GUDHI_CHECK(*vi !=
null_vertex(),
"cannot use the dummy null_vertex() as a real vertex");
677 res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib,
filtration));
678 if (!res_insert.second) {
680 if (res_insert.first->second.filtration() >
filtration) {
682 res_insert.first->second.assign_filtration(
filtration);
686 return std::pair<Simplex_handle, bool>(
null_simplex(),
false);
689 if (static_cast<int>(
simplex.size()) - 1 > dimension_) {
691 dimension_ = static_cast<int>(
simplex.size()) - 1;
720 template<
class InputVertexRange = std::initializer_list<Vertex_handle>>
723 auto first = std::begin(
simplex);
727 return std::pair<Simplex_handle, bool>(
null_simplex(),
true);
730 std::vector<Vertex_handle> copy(first, last);
731 std::sort(std::begin(copy), std::end(copy));
732 return insert_vertex_vector(copy,
filtration);
749 template<
class InputVertexRange = std::initializer_list<Vertex_handle>>
752 auto first = std::begin(Nsimplex);
753 auto last = std::end(Nsimplex);
760 #ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL
762 #endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL
763 std::vector<Vertex_handle> copy;
765 copy.insert(copy.end(), first, last);
766 std::sort(copy.begin(), copy.end());
767 auto last_unique = std::unique(copy.begin(), copy.end());
768 copy.erase(last_unique, copy.end());
771 GUDHI_CHECK(v !=
null_vertex(),
"cannot use the dummy null_vertex() as a real vertex");
774 dimension_ = (std::max)(dimension_, static_cast<int>(std::distance(copy.begin(), copy.end())) - 1);
776 return rec_insert_simplex_and_subfaces_sorted(
root(), copy.begin(), copy.end(),
filtration);
781 template<
class ForwardVertexIterator>
782 std::pair<Simplex_handle, bool> rec_insert_simplex_and_subfaces_sorted(Siblings* sib,
783 ForwardVertexIterator first,
784 ForwardVertexIterator last,
791 auto&& dict = sib->members();
792 auto insertion_result = dict.emplace(vertex_one, Node(sib, filt));
794 bool one_is_new = insertion_result.second;
803 if (++first == last)
return insertion_result;
806 simplex_one->second.assign_children(
new Siblings(sib, vertex_one));
807 auto res = rec_insert_simplex_and_subfaces_sorted(simplex_one->second.children(), first, last, filt);
809 if (res.first !=
null_simplex()) rec_insert_simplex_and_subfaces_sorted(sib, first, last, filt);
817 sh->second.assign_key(
key);
825 return { find_vertex(sh->first), find_vertex(
self_siblings(sh)->parent()) };
829 template<
class SimplexHandle>
831 if (sh->second.children()->parent() == sh->first)
832 return sh->second.children()->oncles();
834 return sh->second.children();
848 dimension_to_be_lowered_ =
false;
863 filtration_vect_.clear();
866 filtration_vect_.push_back(sh);
878 tbb::parallel_sort(filtration_vect_.begin(), filtration_vect_.end(), is_before_in_filtration(
this));
880 std::stable_sort(filtration_vect_.begin(), filtration_vect_.end(), is_before_in_filtration(
this));
897 void rec_coface(std::vector<Vertex_handle> &vertices, Siblings *curr_sib,
int curr_nbVertices,
898 std::vector<Simplex_handle>& cofaces,
bool star,
int nbVertices) {
899 if (!(star || curr_nbVertices <= nbVertices))
902 if (vertices.empty()) {
907 bool addCoface = (star || curr_nbVertices == nbVertices);
911 rec_coface(vertices,
simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
913 if (
simplex->first == vertices.back()) {
915 bool equalDim = (star || curr_nbVertices == nbVertices);
916 bool addCoface = vertices.size() == 1 && equalDim;
923 rec_coface(vertices,
simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
924 vertices.push_back(tmp);
926 }
else if (
simplex->first > vertices.back()) {
931 rec_coface(vertices,
simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
957 assert(codimension >= 0);
959 std::vector<Vertex_handle> copy(rg.begin(), rg.end());
960 if (codimension + static_cast<int>(copy.size()) > dimension_ + 1 ||
961 (codimension == 0 && static_cast<int>(copy.size()) > dimension_))
964 assert(std::is_sorted(copy.begin(), copy.end(), std::greater<Vertex_handle>()));
965 bool star = codimension == 0;
966 rec_coface(copy, &root_, 1, cofaces, star, codimension + static_cast<int>(copy.size()));
983 while (it1 != rg1.end() && it2 != rg2.end()) {
991 return ((it1 == rg1.end()) && (it2 != rg2.end()));
1000 struct is_before_in_filtration {
1006 if (sh1->second.filtration() != sh2->second.filtration()) {
1007 return sh1->second.filtration() < sh2->second.filtration();
1010 return st_->reverse_lexicographic_order(sh1, sh2);
1040 template<
class OneSkeletonGraph>
1045 if (boost::num_vertices(skel_graph) == 0) {
1048 if (num_edges(skel_graph) == 0) {
1054 root_.members_.reserve(boost::num_vertices(skel_graph));
1056 typename boost::graph_traits<OneSkeletonGraph>::vertex_iterator v_it,
1058 for (std::tie(v_it, v_it_end) = boost::vertices(skel_graph); v_it != v_it_end;
1060 root_.members_.emplace_hint(
1061 root_.members_.end(), *v_it,
1062 Node(&root_, boost::get(vertex_filtration_t(), skel_graph, *v_it)));
1064 std::pair<typename boost::graph_traits<OneSkeletonGraph>::edge_iterator,
1065 typename boost::graph_traits<OneSkeletonGraph>::edge_iterator> boost_edges = boost::edges(skel_graph);
1068 for (; boost_edges.first != boost_edges.second; boost_edges.first++) {
1069 auto edge = *(boost_edges.first);
1070 auto u = source(edge, skel_graph);
1071 auto v = target(edge, skel_graph);
1072 if (u == v)
throw "Self-loops are not simplicial";
1080 if (v < u) std::swap(u, v);
1081 auto sh = find_vertex(u);
1083 sh->second.assign_children(
new Siblings(&root_, sh->first));
1086 sh->second.children()->members().emplace(v,
1087 Node(sh->second.children(), boost::get(edge_filtration_t(), skel_graph, edge)));
1103 if (max_dim <= 1)
return;
1104 dimension_ = max_dim;
1105 for (Dictionary_it root_it = root_.members_.begin();
1106 root_it != root_.members_.end(); ++root_it) {
1108 siblings_expansion(root_it->second.children(), max_dim - 1);
1111 dimension_ = max_dim - dimension_;
1116 void siblings_expansion(Siblings * siblings,
1118 if (dimension_ > k) {
1123 Dictionary_it next = siblings->members().begin();
1126 #ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL
1128 #endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL
1129 std::vector<std::pair<Vertex_handle, Node> > inter;
1130 for (Dictionary_it s_h = siblings->members().begin();
1131 s_h != siblings->members().end(); ++s_h, ++next) {
1137 siblings->members().end(),
1138 root_sh->second.children()->members().begin(),
1139 root_sh->second.children()->members().end(),
1140 s_h->second.filtration());
1141 if (inter.size() != 0) {
1142 Siblings * new_sib =
new Siblings(siblings,
1146 s_h->second.assign_children(new_sib);
1147 siblings_expansion(new_sib, k - 1);
1150 s_h->second.assign_children(siblings);
1159 static void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection,
1160 Dictionary_it begin1, Dictionary_it end1,
1161 Dictionary_it begin2, Dictionary_it end2,
1163 if (begin1 == end1 || begin2 == end2)
1166 if (begin1->first == begin2->first) {
1167 Filtration_value filt = (std::max)({begin1->second.filtration(), begin2->second.filtration(), filtration_});
1168 intersection.emplace_back(begin1->first, Node(
nullptr, filt));
1169 if (++begin1 == end1 || ++begin2 == end2)
1171 }
else if (begin1->first < begin2->first) {
1172 if (++begin1 == end1)
1175 if (++begin2 == end2)
1200 template<
typename Blocker >
1203 for (
auto&
simplex : boost::adaptors::reverse(root_.members())) {
1205 siblings_expansion_with_blockers(
simplex.second.children(), max_dim, max_dim - 1, block_simplex);
1212 template<
typename Blocker >
1213 void siblings_expansion_with_blockers(Siblings* siblings,
int max_dim,
int k, Blocker block_simplex) {
1214 if (dimension_ < max_dim - k) {
1215 dimension_ = max_dim - k;
1220 if (siblings->members().size() < 2)
1223 for (
auto simplex = siblings->members().rbegin() + 1;
simplex != siblings->members().rend();
simplex++) {
1224 std::vector<std::pair<Vertex_handle, Node> > intersection;
1225 for(
auto next = siblings->members().rbegin(); next !=
simplex; next++) {
1226 bool to_be_inserted =
true;
1232 to_be_inserted=
false;
1235 filt = (std::max)(filt,
filtration(border_child));
1237 if (to_be_inserted) {
1238 intersection.emplace_back(next->first, Node(
nullptr, filt));
1241 if (intersection.size() != 0) {
1243 Siblings * new_sib =
new Siblings(siblings,
1245 boost::adaptors::reverse(intersection));
1246 std::vector<Vertex_handle> blocked_new_sib_vertex_list;
1248 for (
auto new_sib_member = new_sib->members().begin();
1249 new_sib_member != new_sib->members().end();
1251 bool blocker_result = block_simplex(new_sib_member);
1254 if (blocker_result) {
1255 blocked_new_sib_vertex_list.push_back(new_sib_member->first);
1258 if (blocked_new_sib_vertex_list.size() == new_sib->members().size()) {
1262 simplex->second.assign_children(siblings);
1264 for (
auto& blocked_new_sib_member : blocked_new_sib_vertex_list) {
1265 new_sib->members().erase(blocked_new_sib_member);
1268 simplex->second.assign_children(new_sib);
1269 siblings_expansion_with_blockers(new_sib, max_dim, k - 1, block_simplex);
1273 simplex->second.assign_children(siblings);
1289 if (child == sh->second.children()->members().end())
1307 os <<
key(b_sh) <<
" ";
1322 bool modified =
false;
1324 for (
auto&
simplex : boost::adaptors::reverse(root_.members())) {
1326 modified |= rec_make_filtration_non_decreasing(
simplex.second.children());
1337 bool rec_make_filtration_non_decreasing(Siblings * sib) {
1338 bool modified =
false;
1341 for (
auto&
simplex : boost::adaptors::reverse(sib->members())) {
1350 if (
simplex.second.filtration() < max_filt_border_value) {
1353 simplex.second.assign_filtration(max_filt_border_value);
1356 modified |= rec_make_filtration_non_decreasing(
simplex.second.children());
1380 auto&& list = sib->members();
1381 auto last = std::remove_if(list.begin(), list.end(), [
this,filt](Dit_value_t&
simplex) {
1382 if (
simplex.second.filtration() <= filt)
return false;
1385 dimension_to_be_lowered_ =
true;
1389 bool modified = (last != list.end());
1390 if (last == list.begin() && sib !=
root()) {
1392 sib->oncles()->members()[sib->parent()].assign_children(sib->oncles());
1395 dimension_to_be_lowered_ =
true;
1399 list.erase(last, list.end());
1402 modified |= rec_prune_above_filtration(
simplex.second.children(), filt);
1413 bool lower_upper_bound_dimension() {
1415 dimension_to_be_lowered_ =
false;
1416 int new_dimension = -1;
1421 std::cout <<
" " << vertex;
1423 std::cout << std::endl;
1424 #endif // DEBUG_TRACES
1427 if (sh_dimension >= dimension_)
1430 new_dimension = (std::max)(new_dimension, sh_dimension);
1432 dimension_ = new_dimension;
1450 std::invalid_argument(
"Simplex_tree::remove_maximal_simplex - argument has children"));
1453 Siblings* child = sh->second.children();
1455 if ((child->size() > 1) || (child ==
root())) {
1461 child->oncles()->members().at(child->parent()).assign_children(child->oncles());
1464 dimension_to_be_lowered_ =
true;
1474 std::vector<Simplex_handle> filtration_vect_;
1477 bool dimension_to_be_lowered_ =
false;
1481 template<
typename...T>
1493 template<
typename...T>
1496 std::vector<typename ST::Vertex_handle> simplex;
1497 typename ST::Filtration_value fil;
1502 int dim = static_cast<int> (simplex.size() - 1);
1503 if (max_dim < dim) {
1521 typedef int Vertex_handle;
1522 typedef double Filtration_value;
1523 typedef std::uint32_t Simplex_key;
1524 static const bool store_key =
true;
1525 static const bool store_filtration =
true;
1526 static const bool contiguous_vertices =
false;
1537 typedef int Vertex_handle;
1538 typedef float Filtration_value;
1539 typedef std::uint32_t Simplex_key;
1540 static const bool store_key =
true;
1541 static const bool store_filtration =
true;
1542 static const bool contiguous_vertices =
true;
1549 #endif // SIMPLEX_TREE_H_