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>
27 #include <boost/container/static_vector.hpp>
31 #include <tbb/parallel_sort.h>
40 #include <initializer_list>
61 struct Simplex_tree_options_full_featured;
76 template<
typename SimplexTreeOptions = Simplex_tree_options_full_featured>
95 typedef Simplex_tree_node_explicit_storage<Simplex_tree> Node;
100 typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary;
104 typedef Simplex_tree_siblings<Simplex_tree, Dictionary> Siblings;
108 struct Key_simplex_base_real {
109 Key_simplex_base_real() : key_(-1) {}
115 struct Key_simplex_base_dummy {
116 Key_simplex_base_dummy() {}
121 struct Extended_filtration_data {
124 Extended_filtration_data(){}
127 typedef typename std::conditional<Options::store_key, Key_simplex_base_real, Key_simplex_base_dummy>::type
130 struct Filtration_simplex_base_real {
131 Filtration_simplex_base_real() : filt_(0) {}
137 struct Filtration_simplex_base_dummy {
138 Filtration_simplex_base_dummy() {}
143 Filtration_simplex_base_dummy>::type Filtration_simplex_base;
154 typedef typename Dictionary::iterator Dictionary_it;
155 typedef typename Dictionary_it::value_type Dit_value_t;
157 struct return_first {
221 boost::make_transform_iterator(root_.members_.begin(), return_first()),
222 boost::make_transform_iterator(root_.members_.end(), return_first()));
266 return filtration_vect_;
295 template<
class SimplexHandle>
308 root_(nullptr, null_vertex_),
315 std::clog <<
"Simplex_tree copy constructor" << std::endl;
316 #endif // DEBUG_TRACES
317 copy_from(complex_source);
325 std::clog <<
"Simplex_tree move constructor" << std::endl;
326 #endif // DEBUG_TRACES
327 move_from(complex_source);
331 complex_source.dimension_ = -1;
336 root_members_recursive_deletion();
342 std::clog <<
"Simplex_tree copy assignment" << std::endl;
343 #endif // DEBUG_TRACES
345 if (&complex_source !=
this) {
347 root_members_recursive_deletion();
349 copy_from(complex_source);
359 std::clog <<
"Simplex_tree move assignment" << std::endl;
360 #endif // DEBUG_TRACES
362 if (&complex_source !=
this) {
364 root_members_recursive_deletion();
366 move_from(complex_source);
375 null_vertex_ = complex_source.null_vertex_;
376 filtration_vect_.clear();
377 dimension_ = complex_source.dimension_;
378 auto root_source = complex_source.root_;
381 root_.members() = Dictionary(boost::container::ordered_unique_range, root_source.members().begin(), root_source.members().end());
383 for (
auto& map_el : root_.members()) {
384 map_el.second.assign_children(&root_);
386 rec_copy(&root_, &root_source);
390 void rec_copy(Siblings *sib, Siblings *sib_source) {
391 for (
auto sh = sib->members().begin(), sh_source = sib_source->members().begin();
392 sh != sib->members().end(); ++sh, ++sh_source) {
394 Siblings * newsib =
new Siblings(sib, sh_source->first);
395 newsib->members_.reserve(sh_source->second.children()->members().size());
396 for (
auto & child : sh_source->second.children()->members())
397 newsib->members_.emplace_hint(newsib->members_.end(), child.first, Node(newsib, child.second.filtration()));
398 rec_copy(newsib, sh_source->second.children());
399 sh->second.assign_children(newsib);
406 null_vertex_ = std::move(complex_source.null_vertex_);
407 root_ = std::move(complex_source.root_);
408 filtration_vect_ = std::move(complex_source.filtration_vect_);
409 dimension_ = std::move(complex_source.dimension_);
412 for (
auto& map_el : root_.members()) {
413 if (map_el.second.children() != &(complex_source.root_)) {
415 map_el.second.children()->oncles_ = &root_;
418 GUDHI_CHECK(map_el.second.children()->oncles_ ==
nullptr,
419 std::invalid_argument(
"Simplex_tree move constructor from an invalid Simplex_tree"));
421 map_el.second.assign_children(&root_);
427 void root_members_recursive_deletion() {
428 for (
auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) {
430 rec_delete(sh->second.children());
433 root_.members().clear();
437 void rec_delete(Siblings * sib) {
438 for (
auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
440 rec_delete(sh->second.children());
449 if ((null_vertex_ != st2.null_vertex_) ||
450 (dimension_ != st2.dimension_))
452 return rec_equal(&root_, &st2.root_);
457 return (!(*
this == st2));
462 bool rec_equal(Siblings* s1, Siblings* s2) {
463 if (s1->members().size() != s2->members().size())
465 for (
auto sh1 = s1->members().begin(), sh2 = s2->members().begin();
466 (sh1 != s1->members().end() && sh2 != s2->members().end()); ++sh1, ++sh2) {
467 if (sh1->first != sh2->first || sh1->second.filtration() != sh2->second.filtration())
473 if (!rec_equal(sh1->second.children(), sh2->second.children()))
485 return sh->second.filtration();
495 return sh->second.key();
503 return filtration_vect_[idx];
513 return sh->second.filtration();
515 return std::numeric_limits<Filtration_value>::infinity();
524 std::invalid_argument(
"Simplex_tree::assign_filtration - cannot assign filtration on null_simplex"));
525 sh->second.assign_filtration(fv);
533 return Dictionary_it(
nullptr);
549 return root_.members_.size();
561 auto sib_begin = sib->members().begin();
562 auto sib_end = sib->members().end();
563 size_t simplices_number = sib_end - sib_begin;
564 for (
auto sh = sib_begin; sh != sib_end; ++sh) {
569 return simplices_number;
579 while (curr_sib !=
nullptr) {
581 curr_sib = curr_sib->oncles();
596 if (dimension_to_be_lowered_)
597 lower_upper_bound_dimension();
603 template<
class SimplexHandle>
606 return (sh->second.children()->parent() == sh->first);
616 template<
class InputVertexRange = std::initializer_list<Vertex_handle>>
618 auto first = std::begin(s);
619 auto last = std::end(s);
625 std::vector<Vertex_handle> copy(first, last);
626 std::sort(std::begin(copy), std::end(copy));
627 return find_simplex(copy);
633 Siblings * tmp_sib = &root_;
634 Dictionary_it tmp_dit;
638 GUDHI_CHECK(contiguous_vertices(),
"non-contiguous vertices");
640 if(v < 0 || v >=
static_cast<Vertex_handle>(root_.members_.size()))
642 tmp_dit = root_.members_.begin() + v;
647 tmp_sib = tmp_dit->second.children();
650 tmp_dit = tmp_sib->members_.find(*vi++);
651 if (tmp_dit == tmp_sib->members_.end())
657 tmp_sib = tmp_dit->second.children();
665 assert(contiguous_vertices());
666 return root_.members_.begin() + v;
668 return root_.members_.find(v);
674 bool contiguous_vertices()
const {
675 if (root_.members_.empty())
return true;
676 if (root_.members_.begin()->first != 0)
return false;
677 if (std::prev(root_.members_.end())->first !=
static_cast<Vertex_handle>(root_.members_.size() - 1))
return false;
696 std::pair<Simplex_handle, bool> insert_vertex_vector(
const std::vector<Vertex_handle>&
simplex,
698 Siblings * curr_sib = &root_;
699 std::pair<Simplex_handle, bool> res_insert;
701 for (; vi !=
simplex.end() - 1; ++vi) {
702 GUDHI_CHECK(*vi !=
null_vertex(),
"cannot use the dummy null_vertex() as a real vertex");
703 res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib,
filtration));
705 res_insert.first->second.assign_children(
new Siblings(curr_sib, *vi));
707 curr_sib = res_insert.first->second.children();
709 GUDHI_CHECK(*vi !=
null_vertex(),
"cannot use the dummy null_vertex() as a real vertex");
710 res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib,
filtration));
711 if (!res_insert.second) {
713 if (res_insert.first->second.filtration() >
filtration) {
715 res_insert.first->second.assign_filtration(
filtration);
719 return std::pair<Simplex_handle, bool>(
null_simplex(),
false);
722 if (
static_cast<int>(
simplex.size()) - 1 > dimension_) {
724 dimension_ =
static_cast<int>(
simplex.size()) - 1;
753 template<
class InputVertexRange = std::initializer_list<Vertex_handle>>
756 auto first = std::begin(
simplex);
760 return std::pair<Simplex_handle, bool>(
null_simplex(),
true);
763 std::vector<Vertex_handle> copy(first, last);
764 std::sort(std::begin(copy), std::end(copy));
765 return insert_vertex_vector(copy,
filtration);
782 template<
class InputVertexRange = std::initializer_list<Vertex_handle>>
785 auto first = std::begin(Nsimplex);
786 auto last = std::end(Nsimplex);
791 thread_local std::vector<Vertex_handle> copy;
793 copy.insert(copy.end(), first, last);
794 std::sort(copy.begin(), copy.end());
795 auto last_unique = std::unique(copy.begin(), copy.end());
796 copy.erase(last_unique, copy.end());
799 GUDHI_CHECK(v !=
null_vertex(),
"cannot use the dummy null_vertex() as a real vertex");
802 dimension_ = (std::max)(dimension_,
static_cast<int>(std::distance(copy.begin(), copy.end())) - 1);
804 return rec_insert_simplex_and_subfaces_sorted(
root(), copy.begin(), copy.end(),
filtration);
809 template<
class ForwardVertexIterator>
810 std::pair<Simplex_handle, bool> rec_insert_simplex_and_subfaces_sorted(Siblings* sib,
811 ForwardVertexIterator first,
812 ForwardVertexIterator last,
819 auto&& dict = sib->members();
820 auto insertion_result = dict.emplace(vertex_one, Node(sib, filt));
822 bool one_is_new = insertion_result.second;
831 if (++first == last)
return insertion_result;
834 simplex_one->second.assign_children(
new Siblings(sib, vertex_one));
835 auto res = rec_insert_simplex_and_subfaces_sorted(simplex_one->second.children(), first, last, filt);
837 if (res.first !=
null_simplex()) rec_insert_simplex_and_subfaces_sorted(sib, first, last, filt);
845 sh->second.assign_key(
key);
853 return { find_vertex(sh->first), find_vertex(
self_siblings(sh)->parent()) };
857 template<
class SimplexHandle>
859 if (sh->second.children()->parent() == sh->first)
860 return sh->second.children()->oncles();
862 return sh->second.children();
876 dimension_to_be_lowered_ =
false;
889 filtration_vect_.clear();
892 filtration_vect_.push_back(sh);
904 tbb::parallel_sort(filtration_vect_.begin(), filtration_vect_.end(), is_before_in_filtration(
this));
906 std::stable_sort(filtration_vect_.begin(), filtration_vect_.end(), is_before_in_filtration(
this));
913 if (filtration_vect_.empty()) {
922 filtration_vect_.clear();
938 void rec_coface(std::vector<Vertex_handle> &vertices, Siblings *curr_sib,
int curr_nbVertices,
939 std::vector<Simplex_handle>& cofaces,
bool star,
int nbVertices) {
940 if (!(star || curr_nbVertices <= nbVertices))
943 if (vertices.empty()) {
948 bool addCoface = (star || curr_nbVertices == nbVertices);
952 rec_coface(vertices,
simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
954 if (
simplex->first == vertices.back()) {
956 bool equalDim = (star || curr_nbVertices == nbVertices);
957 bool addCoface = vertices.size() == 1 && equalDim;
964 rec_coface(vertices,
simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
965 vertices.push_back(tmp);
967 }
else if (
simplex->first > vertices.back()) {
972 rec_coface(vertices,
simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
998 assert(codimension >= 0);
1000 std::vector<Vertex_handle> copy(rg.begin(), rg.end());
1001 if (codimension +
static_cast<int>(copy.size()) > dimension_ + 1 ||
1002 (codimension == 0 &&
static_cast<int>(copy.size()) > dimension_))
1005 assert(std::is_sorted(copy.begin(), copy.end(), std::greater<Vertex_handle>()));
1006 bool star = codimension == 0;
1007 rec_coface(copy, &root_, 1, cofaces, star, codimension +
static_cast<int>(copy.size()));
1024 while (it1 != rg1.end() && it2 != rg2.end()) {
1032 return ((it1 == rg1.end()) && (it2 != rg2.end()));
1041 struct is_before_in_filtration {
1047 if (sh1->second.filtration() != sh2->second.filtration()) {
1048 return sh1->second.filtration() < sh2->second.filtration();
1051 return st_->reverse_lexicographic_order(sh1, sh2);
1081 template<
class OneSkeletonGraph>
1086 if (boost::num_vertices(skel_graph) == 0) {
1089 if (num_edges(skel_graph) == 0) {
1095 root_.members_.reserve(boost::num_vertices(skel_graph));
1097 typename boost::graph_traits<OneSkeletonGraph>::vertex_iterator v_it,
1099 for (std::tie(v_it, v_it_end) = boost::vertices(skel_graph); v_it != v_it_end;
1101 root_.members_.emplace_hint(
1102 root_.members_.end(), *v_it,
1103 Node(&root_, boost::get(vertex_filtration_t(), skel_graph, *v_it)));
1105 std::pair<typename boost::graph_traits<OneSkeletonGraph>::edge_iterator,
1106 typename boost::graph_traits<OneSkeletonGraph>::edge_iterator> boost_edges = boost::edges(skel_graph);
1109 for (; boost_edges.first != boost_edges.second; boost_edges.first++) {
1110 auto edge = *(boost_edges.first);
1111 auto u = source(edge, skel_graph);
1112 auto v = target(edge, skel_graph);
1113 if (u == v)
throw "Self-loops are not simplicial";
1121 if (v < u) std::swap(u, v);
1122 auto sh = find_vertex(u);
1124 sh->second.assign_children(
new Siblings(&root_, sh->first));
1127 sh->second.children()->members().emplace(v,
1128 Node(sh->second.children(), boost::get(edge_filtration_t(), skel_graph, edge)));
1144 if (max_dim <= 1)
return;
1146 dimension_ = max_dim;
1147 for (Dictionary_it root_it = root_.members_.begin();
1148 root_it != root_.members_.end(); ++root_it) {
1150 siblings_expansion(root_it->second.children(), max_dim - 1);
1153 dimension_ = max_dim - dimension_;
1158 void siblings_expansion(Siblings * siblings,
1160 if (dimension_ > k) {
1165 Dictionary_it next = siblings->members().begin();
1168 thread_local std::vector<std::pair<Vertex_handle, Node> > inter;
1169 for (Dictionary_it s_h = siblings->members().begin();
1170 s_h != siblings->members().end(); ++s_h, ++next) {
1176 siblings->members().end(),
1177 root_sh->second.children()->members().begin(),
1178 root_sh->second.children()->members().end(),
1179 s_h->second.filtration());
1180 if (inter.size() != 0) {
1181 Siblings * new_sib =
new Siblings(siblings,
1185 s_h->second.assign_children(new_sib);
1186 siblings_expansion(new_sib, k - 1);
1189 s_h->second.assign_children(siblings);
1198 static void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection,
1199 Dictionary_it begin1, Dictionary_it end1,
1200 Dictionary_it begin2, Dictionary_it end2,
1202 if (begin1 == end1 || begin2 == end2)
1205 if (begin1->first == begin2->first) {
1206 Filtration_value filt = (std::max)({begin1->second.filtration(), begin2->second.filtration(), filtration_});
1207 intersection.emplace_back(begin1->first, Node(
nullptr, filt));
1208 if (++begin1 == end1 || ++begin2 == end2)
1210 }
else if (begin1->first < begin2->first) {
1211 if (++begin1 == end1)
1214 if (++begin2 == end2)
1239 template<
typename Blocker >
1242 for (
auto&
simplex : boost::adaptors::reverse(root_.members())) {
1244 siblings_expansion_with_blockers(
simplex.second.children(), max_dim, max_dim - 1, block_simplex);
1251 template<
typename Blocker >
1252 void siblings_expansion_with_blockers(Siblings* siblings,
int max_dim,
int k, Blocker block_simplex) {
1253 if (dimension_ < max_dim - k) {
1254 dimension_ = max_dim - k;
1259 if (siblings->members().size() < 2)
1262 for (
auto simplex = siblings->members().rbegin() + 1;
simplex != siblings->members().rend();
simplex++) {
1263 std::vector<std::pair<Vertex_handle, Node> > intersection;
1264 for(
auto next = siblings->members().rbegin(); next !=
simplex; next++) {
1265 bool to_be_inserted =
true;
1271 to_be_inserted=
false;
1274 filt = (std::max)(filt,
filtration(border_child));
1276 if (to_be_inserted) {
1277 intersection.emplace_back(next->first, Node(
nullptr, filt));
1280 if (intersection.size() != 0) {
1282 Siblings * new_sib =
new Siblings(siblings,
1284 boost::adaptors::reverse(intersection));
1285 std::vector<Vertex_handle> blocked_new_sib_vertex_list;
1287 for (
auto new_sib_member = new_sib->members().begin();
1288 new_sib_member != new_sib->members().end();
1290 bool blocker_result = block_simplex(new_sib_member);
1293 if (blocker_result) {
1294 blocked_new_sib_vertex_list.push_back(new_sib_member->first);
1297 if (blocked_new_sib_vertex_list.size() == new_sib->members().size()) {
1301 simplex->second.assign_children(siblings);
1303 for (
auto& blocked_new_sib_member : blocked_new_sib_vertex_list) {
1304 new_sib->members().erase(blocked_new_sib_member);
1307 simplex->second.assign_children(new_sib);
1308 siblings_expansion_with_blockers(new_sib, max_dim, k - 1, block_simplex);
1312 simplex->second.assign_children(siblings);
1328 if (child == sh->second.children()->members().end())
1346 os <<
key(b_sh) <<
" ";
1360 bool modified =
false;
1362 for (
auto&
simplex : boost::adaptors::reverse(root_.members())) {
1364 modified |= rec_make_filtration_non_decreasing(
simplex.second.children());
1377 bool rec_make_filtration_non_decreasing(Siblings * sib) {
1378 bool modified =
false;
1381 for (
auto&
simplex : boost::adaptors::reverse(sib->members())) {
1392 if (!(
simplex.second.filtration() >= max_filt_border_value)) {
1395 simplex.second.assign_filtration(max_filt_border_value);
1398 modified |= rec_make_filtration_non_decreasing(
simplex.second.children());
1422 auto&& list = sib->members();
1423 auto last = std::remove_if(list.begin(), list.end(), [
this,filt](Dit_value_t&
simplex) {
1424 if (simplex.second.filtration() <= filt) return false;
1425 if (has_children(&simplex)) rec_delete(simplex.second.children());
1427 dimension_to_be_lowered_ = true;
1431 bool modified = (last != list.end());
1432 if (last == list.begin() && sib !=
root()) {
1434 sib->oncles()->members()[sib->parent()].assign_children(sib->oncles());
1437 dimension_to_be_lowered_ =
true;
1441 list.erase(last, list.end());
1444 modified |= rec_prune_above_filtration(
simplex.second.children(), filt);
1455 bool lower_upper_bound_dimension() {
1457 dimension_to_be_lowered_ =
false;
1458 int new_dimension = -1;
1463 std::clog <<
" " << vertex;
1465 std::clog << std::endl;
1466 #endif // DEBUG_TRACES
1469 if (sh_dimension >= dimension_)
1472 new_dimension = (std::max)(new_dimension, sh_dimension);
1474 dimension_ = new_dimension;
1491 std::invalid_argument(
"Simplex_tree::remove_maximal_simplex - argument has children"));
1494 Siblings* child = sh->second.children();
1496 if ((child->size() > 1) || (child ==
root())) {
1502 child->oncles()->members().at(child->parent()).assign_children(child->oncles());
1505 dimension_to_be_lowered_ =
true;
1526 std::pair<Filtration_value, Extended_simplex_type> p;
1529 if (f >= -2 && f <= -1){
1530 p.first = minval + (maxval-minval)*(f + 2); p.second = Extended_simplex_type::UP;
1532 else if (f >= 1 && f <= 2){
1533 p.first = minval - (maxval-minval)*(f - 2); p.second = Extended_simplex_type::DOWN;
1536 p.first = std::numeric_limits<Filtration_value>::quiet_NaN(); p.second = Extended_simplex_type::EXTRA;
1558 Vertex_handle maxvert = std::numeric_limits<Vertex_handle>::min();
1559 Filtration_value minval = std::numeric_limits<Filtration_value>::infinity();
1560 Filtration_value maxval = -std::numeric_limits<Filtration_value>::infinity();
1561 for (
auto sh = root_.members().begin(); sh != root_.members().end(); ++sh){
1563 minval = std::min(minval, f);
1564 maxval = std::max(maxval, f);
1565 maxvert = std::max(sh->first, maxvert);
1568 GUDHI_CHECK(maxvert < std::numeric_limits<Vertex_handle>::max(), std::invalid_argument(
"Simplex_tree contains a vertex with the largest Vertex_handle"));
1577 std::vector<Vertex_handle> vr;
1585 auto sh = this->
find(vr);
1588 vr.push_back(maxvert);
1608 Extended_filtration_data efd(minval, maxval);
1617 auto filt = filtration_(sh);
1619 if(filtration_(find_vertex(v)) == filt)
1633 auto end = std::end(vertices);
1634 auto vi = std::begin(vertices);
1635 GUDHI_CHECK(vi != end,
"empty simplex");
1638 GUDHI_CHECK(vi != end,
"simplex of dimension 0");
1639 if(std::next(vi) == end)
return sh;
1640 boost::container::static_vector<Vertex_handle, 40> suffix;
1641 suffix.push_back(v0);
1642 auto filt = filtration_(sh);
1646 auto&& children1 = find_vertex(v)->second.children()->members_;
1647 for(
auto w : suffix){
1650 if(filtration_(s) == filt)
1653 suffix.push_back(v);
1664 auto filt = filtration_(sh);
1667 if(filtration_(b) == filt)
1678 std::vector<Simplex_handle> filtration_vect_;
1681 bool dimension_to_be_lowered_ =
false;
1685 template<
typename...T>
1697 template<
typename...T>
1700 std::vector<typename ST::Vertex_handle> simplex;
1701 typename ST::Filtration_value fil;
1706 int dim =
static_cast<int> (simplex.size() - 1);
1707 if (max_dim < dim) {
1725 typedef int Vertex_handle;
1726 typedef double Filtration_value;
1727 typedef std::uint32_t Simplex_key;
1728 static const bool store_key =
true;
1729 static const bool store_filtration =
true;
1730 static const bool contiguous_vertices =
false;
1741 typedef int Vertex_handle;
1742 typedef float Filtration_value;
1743 typedef std::uint32_t Simplex_key;
1744 static const bool store_key =
true;
1745 static const bool store_filtration =
true;
1746 static const bool contiguous_vertices =
true;
1753 #endif // SIMPLEX_TREE_H_