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>
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>
27 #include <boost/container/static_vector.hpp>
30 #include <tbb/parallel_sort.h>
38 #include <initializer_list>
59 struct Simplex_tree_options_full_featured;
74 template<
typename SimplexTreeOptions = Simplex_tree_options_full_featured>
93 typedef Simplex_tree_node_explicit_storage<Simplex_tree> Node;
98 typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary;
102 typedef Simplex_tree_siblings<Simplex_tree, Dictionary> Siblings;
106 struct Key_simplex_base_real {
107 Key_simplex_base_real() : key_(-1) {}
113 struct Key_simplex_base_dummy {
114 Key_simplex_base_dummy() {}
119 struct Extended_filtration_data {
122 Extended_filtration_data(){}
125 typedef typename std::conditional<Options::store_key, Key_simplex_base_real, Key_simplex_base_dummy>::type
128 struct Filtration_simplex_base_real {
129 Filtration_simplex_base_real() : filt_(0) {}
135 struct Filtration_simplex_base_dummy {
136 Filtration_simplex_base_dummy() {}
141 Filtration_simplex_base_dummy>::type Filtration_simplex_base;
152 typedef typename Dictionary::iterator Dictionary_it;
153 typedef typename Dictionary_it::value_type Dit_value_t;
155 struct return_first {
219 boost::make_transform_iterator(root_.members_.begin(), return_first()),
220 boost::make_transform_iterator(root_.members_.end(), return_first()));
264 return filtration_vect_;
293 template<
class SimplexHandle>
306 root_(nullptr, null_vertex_),
313 std::clog <<
"Simplex_tree copy constructor" << std::endl;
314 #endif // DEBUG_TRACES
315 copy_from(complex_source);
323 std::clog <<
"Simplex_tree move constructor" << std::endl;
324 #endif // DEBUG_TRACES
325 move_from(complex_source);
329 complex_source.dimension_ = -1;
334 root_members_recursive_deletion();
340 std::clog <<
"Simplex_tree copy assignment" << std::endl;
341 #endif // DEBUG_TRACES
343 if (&complex_source !=
this) {
345 root_members_recursive_deletion();
347 copy_from(complex_source);
357 std::clog <<
"Simplex_tree move assignment" << std::endl;
358 #endif // DEBUG_TRACES
360 if (&complex_source !=
this) {
362 root_members_recursive_deletion();
364 move_from(complex_source);
373 null_vertex_ = complex_source.null_vertex_;
374 filtration_vect_.clear();
375 dimension_ = complex_source.dimension_;
376 auto root_source = complex_source.root_;
379 root_.members() = Dictionary(boost::container::ordered_unique_range, root_source.members().begin(), root_source.members().end());
381 for (
auto& map_el : root_.members()) {
382 map_el.second.assign_children(&root_);
384 rec_copy(&root_, &root_source);
388 void rec_copy(Siblings *sib, Siblings *sib_source) {
389 for (
auto sh = sib->members().begin(), sh_source = sib_source->members().begin();
390 sh != sib->members().end(); ++sh, ++sh_source) {
392 Siblings * newsib =
new Siblings(sib, sh_source->first);
393 newsib->members_.reserve(sh_source->second.children()->members().size());
394 for (
auto & child : sh_source->second.children()->members())
395 newsib->members_.emplace_hint(newsib->members_.end(), child.first, Node(newsib, child.second.filtration()));
396 rec_copy(newsib, sh_source->second.children());
397 sh->second.assign_children(newsib);
404 null_vertex_ = std::move(complex_source.null_vertex_);
405 root_ = std::move(complex_source.root_);
406 filtration_vect_ = std::move(complex_source.filtration_vect_);
407 dimension_ = std::move(complex_source.dimension_);
410 for (
auto& map_el : root_.members()) {
411 if (map_el.second.children() != &(complex_source.root_)) {
413 map_el.second.children()->oncles_ = &root_;
416 GUDHI_CHECK(map_el.second.children()->oncles_ ==
nullptr,
417 std::invalid_argument(
"Simplex_tree move constructor from an invalid Simplex_tree"));
419 map_el.second.assign_children(&root_);
425 void root_members_recursive_deletion() {
426 for (
auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) {
428 rec_delete(sh->second.children());
431 root_.members().clear();
435 void rec_delete(Siblings * sib) {
436 for (
auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
438 rec_delete(sh->second.children());
447 if ((null_vertex_ != st2.null_vertex_) ||
448 (dimension_ != st2.dimension_))
450 return rec_equal(&root_, &st2.root_);
455 return (!(*
this == st2));
460 bool rec_equal(Siblings* s1, Siblings* s2) {
461 if (s1->members().size() != s2->members().size())
463 for (
auto sh1 = s1->members().begin(), sh2 = s2->members().begin();
464 (sh1 != s1->members().end() && sh2 != s2->members().end()); ++sh1, ++sh2) {
465 if (sh1->first != sh2->first || sh1->second.filtration() != sh2->second.filtration())
471 if (!rec_equal(sh1->second.children(), sh2->second.children()))
483 return sh->second.filtration();
493 return sh->second.key();
501 return filtration_vect_[idx];
511 return sh->second.filtration();
513 return std::numeric_limits<Filtration_value>::infinity();
522 std::invalid_argument(
"Simplex_tree::assign_filtration - cannot assign filtration on null_simplex"));
523 sh->second.assign_filtration(fv);
531 return Dictionary_it(
nullptr);
547 return root_.members_.size();
559 auto sib_begin = sib->members().begin();
560 auto sib_end = sib->members().end();
561 size_t simplices_number = sib_end - sib_begin;
562 for (
auto sh = sib_begin; sh != sib_end; ++sh) {
567 return simplices_number;
577 while (curr_sib !=
nullptr) {
579 curr_sib = curr_sib->oncles();
594 if (dimension_to_be_lowered_)
595 lower_upper_bound_dimension();
601 template<
class SimplexHandle>
604 return (sh->second.children()->parent() == sh->first);
614 template<
class InputVertexRange = std::initializer_list<Vertex_handle>>
616 auto first = std::begin(s);
617 auto last = std::end(s);
623 std::vector<Vertex_handle> copy(first, last);
624 std::sort(std::begin(copy), std::end(copy));
625 return find_simplex(copy);
631 Siblings * tmp_sib = &root_;
632 Dictionary_it tmp_dit;
636 GUDHI_CHECK(contiguous_vertices(),
"non-contiguous vertices");
638 if(v < 0 || v >=
static_cast<Vertex_handle>(root_.members_.size()))
640 tmp_dit = root_.members_.begin() + v;
645 tmp_sib = tmp_dit->second.children();
648 tmp_dit = tmp_sib->members_.find(*vi++);
649 if (tmp_dit == tmp_sib->members_.end())
655 tmp_sib = tmp_dit->second.children();
663 assert(contiguous_vertices());
664 return root_.members_.begin() + v;
666 return root_.members_.find(v);
672 bool contiguous_vertices()
const {
673 if (root_.members_.empty())
return true;
674 if (root_.members_.begin()->first != 0)
return false;
675 if (std::prev(root_.members_.end())->first !=
static_cast<Vertex_handle>(root_.members_.size() - 1))
return false;
694 std::pair<Simplex_handle, bool> insert_vertex_vector(
const std::vector<Vertex_handle>&
simplex,
696 Siblings * curr_sib = &root_;
697 std::pair<Simplex_handle, bool> res_insert;
699 for (; vi !=
simplex.end() - 1; ++vi) {
700 GUDHI_CHECK(*vi !=
null_vertex(),
"cannot use the dummy null_vertex() as a real vertex");
701 res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib,
filtration));
703 res_insert.first->second.assign_children(
new Siblings(curr_sib, *vi));
705 curr_sib = res_insert.first->second.children();
707 GUDHI_CHECK(*vi !=
null_vertex(),
"cannot use the dummy null_vertex() as a real vertex");
708 res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib,
filtration));
709 if (!res_insert.second) {
711 if (res_insert.first->second.filtration() >
filtration) {
713 res_insert.first->second.assign_filtration(
filtration);
717 return std::pair<Simplex_handle, bool>(
null_simplex(),
false);
720 if (
static_cast<int>(
simplex.size()) - 1 > dimension_) {
722 dimension_ =
static_cast<int>(
simplex.size()) - 1;
751 template<
class InputVertexRange = std::initializer_list<Vertex_handle>>
754 auto first = std::begin(
simplex);
758 return std::pair<Simplex_handle, bool>(
null_simplex(),
true);
761 std::vector<Vertex_handle> copy(first, last);
762 std::sort(std::begin(copy), std::end(copy));
763 return insert_vertex_vector(copy,
filtration);
780 template<
class InputVertexRange = std::initializer_list<Vertex_handle>>
783 auto first = std::begin(Nsimplex);
784 auto last = std::end(Nsimplex);
789 thread_local std::vector<Vertex_handle> copy;
791 copy.insert(copy.end(), first, last);
792 std::sort(copy.begin(), copy.end());
793 auto last_unique = std::unique(copy.begin(), copy.end());
794 copy.erase(last_unique, copy.end());
797 GUDHI_CHECK(v !=
null_vertex(),
"cannot use the dummy null_vertex() as a real vertex");
800 dimension_ = (std::max)(dimension_,
static_cast<int>(std::distance(copy.begin(), copy.end())) - 1);
802 return rec_insert_simplex_and_subfaces_sorted(
root(), copy.begin(), copy.end(),
filtration);
807 template<
class ForwardVertexIterator>
808 std::pair<Simplex_handle, bool> rec_insert_simplex_and_subfaces_sorted(Siblings* sib,
809 ForwardVertexIterator first,
810 ForwardVertexIterator last,
817 auto&& dict = sib->members();
818 auto insertion_result = dict.emplace(vertex_one, Node(sib, filt));
820 bool one_is_new = insertion_result.second;
829 if (++first == last)
return insertion_result;
832 simplex_one->second.assign_children(
new Siblings(sib, vertex_one));
833 auto res = rec_insert_simplex_and_subfaces_sorted(simplex_one->second.children(), first, last, filt);
835 if (res.first !=
null_simplex()) rec_insert_simplex_and_subfaces_sorted(sib, first, last, filt);
843 sh->second.assign_key(
key);
851 return { find_vertex(sh->first), find_vertex(
self_siblings(sh)->parent()) };
855 template<
class SimplexHandle>
857 if (sh->second.children()->parent() == sh->first)
858 return sh->second.children()->oncles();
860 return sh->second.children();
874 dimension_to_be_lowered_ =
false;
887 filtration_vect_.clear();
890 filtration_vect_.push_back(sh);
902 tbb::parallel_sort(filtration_vect_.begin(), filtration_vect_.end(), is_before_in_filtration(
this));
904 std::stable_sort(filtration_vect_.begin(), filtration_vect_.end(), is_before_in_filtration(
this));
911 if (filtration_vect_.empty()) {
920 filtration_vect_.clear();
936 void rec_coface(std::vector<Vertex_handle> &vertices, Siblings *curr_sib,
int curr_nbVertices,
937 std::vector<Simplex_handle>& cofaces,
bool star,
int nbVertices) {
938 if (!(star || curr_nbVertices <= nbVertices))
941 if (vertices.empty()) {
946 bool addCoface = (star || curr_nbVertices == nbVertices);
950 rec_coface(vertices,
simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
952 if (
simplex->first == vertices.back()) {
954 bool equalDim = (star || curr_nbVertices == nbVertices);
955 bool addCoface = vertices.size() == 1 && equalDim;
962 rec_coface(vertices,
simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
963 vertices.push_back(tmp);
965 }
else if (
simplex->first > vertices.back()) {
970 rec_coface(vertices,
simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
996 assert(codimension >= 0);
998 std::vector<Vertex_handle> copy(rg.begin(), rg.end());
999 if (codimension +
static_cast<int>(copy.size()) > dimension_ + 1 ||
1000 (codimension == 0 &&
static_cast<int>(copy.size()) > dimension_))
1003 assert(std::is_sorted(copy.begin(), copy.end(), std::greater<Vertex_handle>()));
1004 bool star = codimension == 0;
1005 rec_coface(copy, &root_, 1, cofaces, star, codimension +
static_cast<int>(copy.size()));
1022 while (it1 != rg1.end() && it2 != rg2.end()) {
1030 return ((it1 == rg1.end()) && (it2 != rg2.end()));
1039 struct is_before_in_filtration {
1045 if (sh1->second.filtration() != sh2->second.filtration()) {
1046 return sh1->second.filtration() < sh2->second.filtration();
1049 return st_->reverse_lexicographic_order(sh1, sh2);
1079 template<
class OneSkeletonGraph>
1084 if (boost::num_vertices(skel_graph) == 0) {
1087 if (num_edges(skel_graph) == 0) {
1093 root_.members_.reserve(boost::num_vertices(skel_graph));
1095 typename boost::graph_traits<OneSkeletonGraph>::vertex_iterator v_it,
1097 for (std::tie(v_it, v_it_end) = boost::vertices(skel_graph); v_it != v_it_end;
1099 root_.members_.emplace_hint(
1100 root_.members_.end(), *v_it,
1101 Node(&root_, boost::get(vertex_filtration_t(), skel_graph, *v_it)));
1103 std::pair<typename boost::graph_traits<OneSkeletonGraph>::edge_iterator,
1104 typename boost::graph_traits<OneSkeletonGraph>::edge_iterator> boost_edges = boost::edges(skel_graph);
1107 for (; boost_edges.first != boost_edges.second; boost_edges.first++) {
1108 auto edge = *(boost_edges.first);
1109 auto u = source(edge, skel_graph);
1110 auto v = target(edge, skel_graph);
1111 if (u == v)
throw "Self-loops are not simplicial";
1119 if (v < u) std::swap(u, v);
1120 auto sh = find_vertex(u);
1122 sh->second.assign_children(
new Siblings(&root_, sh->first));
1125 sh->second.children()->members().emplace(v,
1126 Node(sh->second.children(), boost::get(edge_filtration_t(), skel_graph, edge)));
1142 if (max_dim <= 1)
return;
1144 dimension_ = max_dim;
1145 for (Dictionary_it root_it = root_.members_.begin();
1146 root_it != root_.members_.end(); ++root_it) {
1148 siblings_expansion(root_it->second.children(), max_dim - 1);
1151 dimension_ = max_dim - dimension_;
1156 void siblings_expansion(Siblings * siblings,
1158 if (dimension_ > k) {
1163 Dictionary_it next = siblings->members().begin();
1166 thread_local std::vector<std::pair<Vertex_handle, Node> > inter;
1167 for (Dictionary_it s_h = siblings->members().begin();
1168 s_h != siblings->members().end(); ++s_h, ++next) {
1174 siblings->members().end(),
1175 root_sh->second.children()->members().begin(),
1176 root_sh->second.children()->members().end(),
1177 s_h->second.filtration());
1178 if (inter.size() != 0) {
1179 Siblings * new_sib =
new Siblings(siblings,
1183 s_h->second.assign_children(new_sib);
1184 siblings_expansion(new_sib, k - 1);
1187 s_h->second.assign_children(siblings);
1196 static void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection,
1197 Dictionary_it begin1, Dictionary_it end1,
1198 Dictionary_it begin2, Dictionary_it end2,
1200 if (begin1 == end1 || begin2 == end2)
1203 if (begin1->first == begin2->first) {
1204 Filtration_value filt = (std::max)({begin1->second.filtration(), begin2->second.filtration(), filtration_});
1205 intersection.emplace_back(begin1->first, Node(
nullptr, filt));
1206 if (++begin1 == end1 || ++begin2 == end2)
1208 }
else if (begin1->first < begin2->first) {
1209 if (++begin1 == end1)
1212 if (++begin2 == end2)
1237 template<
typename Blocker >
1240 for (
auto&
simplex : boost::adaptors::reverse(root_.members())) {
1242 siblings_expansion_with_blockers(
simplex.second.children(), max_dim, max_dim - 1, block_simplex);
1249 template<
typename Blocker >
1250 void siblings_expansion_with_blockers(Siblings* siblings,
int max_dim,
int k, Blocker block_simplex) {
1251 if (dimension_ < max_dim - k) {
1252 dimension_ = max_dim - k;
1257 if (siblings->members().size() < 2)
1260 for (
auto simplex = siblings->members().rbegin() + 1;
simplex != siblings->members().rend();
simplex++) {
1261 std::vector<std::pair<Vertex_handle, Node> > intersection;
1262 for(
auto next = siblings->members().rbegin(); next !=
simplex; next++) {
1263 bool to_be_inserted =
true;
1269 to_be_inserted=
false;
1272 filt = (std::max)(filt,
filtration(border_child));
1274 if (to_be_inserted) {
1275 intersection.emplace_back(next->first, Node(
nullptr, filt));
1278 if (intersection.size() != 0) {
1280 Siblings * new_sib =
new Siblings(siblings,
1282 boost::adaptors::reverse(intersection));
1283 std::vector<Vertex_handle> blocked_new_sib_vertex_list;
1285 for (
auto new_sib_member = new_sib->members().begin();
1286 new_sib_member != new_sib->members().end();
1288 bool blocker_result = block_simplex(new_sib_member);
1291 if (blocker_result) {
1292 blocked_new_sib_vertex_list.push_back(new_sib_member->first);
1295 if (blocked_new_sib_vertex_list.size() == new_sib->members().size()) {
1299 simplex->second.assign_children(siblings);
1301 for (
auto& blocked_new_sib_member : blocked_new_sib_vertex_list) {
1302 new_sib->members().erase(blocked_new_sib_member);
1305 simplex->second.assign_children(new_sib);
1306 siblings_expansion_with_blockers(new_sib, max_dim, k - 1, block_simplex);
1310 simplex->second.assign_children(siblings);
1326 if (child == sh->second.children()->members().end())
1344 os <<
key(b_sh) <<
" ";
1358 bool modified =
false;
1360 for (
auto&
simplex : boost::adaptors::reverse(root_.members())) {
1362 modified |= rec_make_filtration_non_decreasing(
simplex.second.children());
1375 bool rec_make_filtration_non_decreasing(Siblings * sib) {
1376 bool modified =
false;
1379 for (
auto&
simplex : boost::adaptors::reverse(sib->members())) {
1390 if (!(
simplex.second.filtration() >= max_filt_border_value)) {
1393 simplex.second.assign_filtration(max_filt_border_value);
1396 modified |= rec_make_filtration_non_decreasing(
simplex.second.children());
1420 auto&& list = sib->members();
1421 auto last = std::remove_if(list.begin(), list.end(), [
this,filt](Dit_value_t&
simplex) {
1422 if (simplex.second.filtration() <= filt) return false;
1423 if (has_children(&simplex)) rec_delete(simplex.second.children());
1425 dimension_to_be_lowered_ = true;
1429 bool modified = (last != list.end());
1430 if (last == list.begin() && sib !=
root()) {
1432 sib->oncles()->members()[sib->parent()].assign_children(sib->oncles());
1435 dimension_to_be_lowered_ =
true;
1439 list.erase(last, list.end());
1442 modified |= rec_prune_above_filtration(
simplex.second.children(), filt);
1453 bool lower_upper_bound_dimension() {
1455 dimension_to_be_lowered_ =
false;
1456 int new_dimension = -1;
1461 std::clog <<
" " << vertex;
1463 std::clog << std::endl;
1464 #endif // DEBUG_TRACES
1467 if (sh_dimension >= dimension_)
1470 new_dimension = (std::max)(new_dimension, sh_dimension);
1472 dimension_ = new_dimension;
1489 std::invalid_argument(
"Simplex_tree::remove_maximal_simplex - argument has children"));
1492 Siblings* child = sh->second.children();
1494 if ((child->size() > 1) || (child ==
root())) {
1500 child->oncles()->members().at(child->parent()).assign_children(child->oncles());
1503 dimension_to_be_lowered_ =
true;
1524 std::pair<Filtration_value, Extended_simplex_type> p;
1527 if (f >= -2 && f <= -1){
1528 p.first = minval + (maxval-minval)*(f + 2); p.second = Extended_simplex_type::UP;
1530 else if (f >= 1 && f <= 2){
1531 p.first = minval - (maxval-minval)*(f - 2); p.second = Extended_simplex_type::DOWN;
1534 p.first = std::numeric_limits<Filtration_value>::quiet_NaN(); p.second = Extended_simplex_type::EXTRA;
1556 Vertex_handle maxvert = std::numeric_limits<Vertex_handle>::min();
1557 Filtration_value minval = std::numeric_limits<Filtration_value>::infinity();
1558 Filtration_value maxval = -std::numeric_limits<Filtration_value>::infinity();
1559 for (
auto sh = root_.members().begin(); sh != root_.members().end(); ++sh){
1561 minval = std::min(minval, f);
1562 maxval = std::max(maxval, f);
1563 maxvert = std::max(sh->first, maxvert);
1566 GUDHI_CHECK(maxvert < std::numeric_limits<Vertex_handle>::max(), std::invalid_argument(
"Simplex_tree contains a vertex with the largest Vertex_handle"));
1575 std::vector<Vertex_handle> vr;
1583 auto sh = this->
find(vr);
1586 vr.push_back(maxvert);
1606 Extended_filtration_data efd(minval, maxval);
1615 auto filt = filtration_(sh);
1617 if(filtration_(find_vertex(v)) == filt)
1631 auto end = std::end(vertices);
1632 auto vi = std::begin(vertices);
1633 GUDHI_CHECK(vi != end,
"empty simplex");
1636 GUDHI_CHECK(vi != end,
"simplex of dimension 0");
1637 if(std::next(vi) == end)
return sh;
1638 boost::container::static_vector<Vertex_handle, 40> suffix;
1639 suffix.push_back(v0);
1640 auto filt = filtration_(sh);
1644 auto&& children1 = find_vertex(v)->second.children()->members_;
1645 for(
auto w : suffix){
1648 if(filtration_(s) == filt)
1651 suffix.push_back(v);
1662 auto filt = filtration_(sh);
1665 if(filtration_(b) == filt)
1680 rec_reset_filtration(&root_, filt_value, min_dim);
1690 void rec_reset_filtration(Siblings * sib,
Filtration_value filt_value,
int min_depth) {
1691 for (
auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
1692 if (min_depth <= 0) {
1693 sh->second.assign_filtration(filt_value);
1696 rec_reset_filtration(sh->second.children(), filt_value, min_depth - 1);
1707 std::vector<Simplex_handle> filtration_vect_;
1710 bool dimension_to_be_lowered_ =
false;
1714 template<
typename...T>
1726 template<
typename...T>
1729 std::vector<typename ST::Vertex_handle> simplex;
1730 typename ST::Filtration_value fil;
1735 int dim =
static_cast<int> (simplex.size() - 1);
1736 if (max_dim < dim) {
1754 typedef int Vertex_handle;
1755 typedef double Filtration_value;
1756 typedef std::uint32_t Simplex_key;
1757 static const bool store_key =
true;
1758 static const bool store_filtration =
true;
1759 static const bool contiguous_vertices =
false;
1770 typedef int Vertex_handle;
1771 typedef float Filtration_value;
1772 typedef std::uint32_t Simplex_key;
1773 static const bool store_key =
true;
1774 static const bool store_filtration =
true;
1775 static const bool contiguous_vertices =
true;
1782 #endif // SIMPLEX_TREE_H_