Loading...
Searching...
No Matches
Simplex_tree.h
1/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2 * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3 * Author(s): Clément Maria
4 *
5 * Copyright (C) 2014 Inria
6 *
7 * Modification(s):
8 * - Vincent Rouvreau: Add de/serialize methods for pickle feature
9 * - YYYY/MM Author: Description of the modification
10 */
11
12#ifndef SIMPLEX_TREE_H_
13#define SIMPLEX_TREE_H_
14
15#include <gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h>
16#include <gudhi/Simplex_tree/Simplex_tree_siblings.h>
17#include <gudhi/Simplex_tree/Simplex_tree_iterators.h>
18#include <gudhi/Simplex_tree/indexing_tag.h>
19#include <gudhi/Simplex_tree/serialization_utils.h> // for Gudhi::simplex_tree::de/serialize_trivial
20
21#include <gudhi/reader_utils.h>
23#include <gudhi/Debug_utils.h>
24
25#include <boost/container/flat_map.hpp>
26#include <boost/iterator/transform_iterator.hpp>
27#include <boost/graph/adjacency_list.hpp>
28#include <boost/range/adaptor/reversed.hpp>
29#include <boost/range/adaptor/transformed.hpp>
30#include <boost/range/size.hpp>
31#include <boost/container/static_vector.hpp>
32
33#ifdef GUDHI_USE_TBB
34#include <tbb/parallel_sort.h>
35#endif
36
37#include <utility>
38#include <vector>
39#include <functional> // for greater<>
40#include <stdexcept>
41#include <limits> // Inf
42#include <initializer_list>
43#include <algorithm> // for std::max
44#include <cstdint> // for std::uint32_t
45#include <iterator> // for std::distance
46
47namespace Gudhi {
48
65enum class Extended_simplex_type {UP, DOWN, EXTRA};
66
67struct Simplex_tree_options_full_featured;
68
82template<typename SimplexTreeOptions = Simplex_tree_options_full_featured>
84 public:
86 typedef typename Options::Indexing_tag Indexing_tag;
99
100 /* Type of node in the simplex tree. */
102 /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */
103 // Note: this wastes space when Vertex_handle is 32 bits and Node is aligned on 64 bits. It would be better to use a
104 // flat_set (with our own comparator) where we can control the layout of the struct (put Vertex_handle and
105 // Simplex_key next to each other).
106 typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary;
107
110
111
112
113 struct Key_simplex_base_real {
114 Key_simplex_base_real() : key_(-1) {}
115 void assign_key(Simplex_key k) { key_ = k; }
116 Simplex_key key() const { return key_; }
117 private:
118 Simplex_key key_;
119 };
120 struct Key_simplex_base_dummy {
121 Key_simplex_base_dummy() {}
122 // Undefined so it will not link
123 void assign_key(Simplex_key);
124 Simplex_key key() const;
125 };
126 struct Extended_filtration_data {
127 Filtration_value minval;
128 Filtration_value maxval;
129 Extended_filtration_data(){}
130 Extended_filtration_data(Filtration_value vmin, Filtration_value vmax): minval(vmin), maxval(vmax) {}
131 };
132 typedef typename std::conditional<Options::store_key, Key_simplex_base_real, Key_simplex_base_dummy>::type
133 Key_simplex_base;
134
135 struct Filtration_simplex_base_real {
136 Filtration_simplex_base_real() : filt_(0) {}
137 void assign_filtration(Filtration_value f) { filt_ = f; }
138 Filtration_value filtration() const { return filt_; }
139 private:
140 Filtration_value filt_;
141 };
142 struct Filtration_simplex_base_dummy {
143 Filtration_simplex_base_dummy() {}
144 void assign_filtration(Filtration_value GUDHI_CHECK_code(f)) { GUDHI_CHECK(f == 0, "filtration value specified for a complex that does not store them"); }
145 Filtration_value filtration() const { return 0; }
146 };
147 typedef typename std::conditional<Options::store_filtration, Filtration_simplex_base_real,
148 Filtration_simplex_base_dummy>::type Filtration_simplex_base;
149
150 public:
156 typedef typename Dictionary::iterator Simplex_handle;
157
158 private:
159 typedef typename Dictionary::iterator Dictionary_it;
160 typedef typename Dictionary_it::value_type Dit_value_t;
161
162 struct return_first {
163 Vertex_handle operator()(const Dit_value_t& p_sh) const {
164 return p_sh.first;
165 }
166 };
167
168 public:
180 typedef boost::transform_iterator<return_first, Dictionary_it> Complex_vertex_iterator;
182 typedef boost::iterator_range<Complex_vertex_iterator> Complex_vertex_range;
188 typedef boost::iterator_range<Simplex_vertex_iterator> Simplex_vertex_range;
190 typedef std::vector<Simplex_handle> Cofaces_simplex_range;
196 typedef boost::iterator_range<Boundary_simplex_iterator> Boundary_simplex_range;
202 typedef boost::iterator_range<Boundary_opposite_vertex_simplex_iterator> Boundary_opposite_vertex_simplex_range;
208 typedef boost::iterator_range<Complex_simplex_iterator> Complex_simplex_range;
216 typedef boost::iterator_range<Skeleton_simplex_iterator> Skeleton_simplex_range;
218 typedef std::vector<Simplex_handle> Filtration_simplex_range;
222 typedef typename Filtration_simplex_range::const_iterator Filtration_simplex_iterator;
223
224 /* @} */ // end name range and iterator types
232 boost::make_transform_iterator(root_.members_.begin(), return_first()),
233 boost::make_transform_iterator(root_.members_.end(), return_first()));
234 }
235
244 }
245
258 }
259
277 return filtration_vect_;
278 }
279
287 GUDHI_CHECK(sh != null_simplex(), "empty simplex");
290 }
291
306 template<class SimplexHandle>
310 }
311
323 template<class SimplexHandle>
327 }
328 // end range and iterator methods
335 : null_vertex_(-1),
336 root_(nullptr, null_vertex_),
337 filtration_vect_(),
338 dimension_(-1) { }
339
341 Simplex_tree(const Simplex_tree& complex_source) {
342#ifdef DEBUG_TRACES
343 std::clog << "Simplex_tree copy constructor" << std::endl;
344#endif // DEBUG_TRACES
345 copy_from(complex_source);
346 }
347
351 Simplex_tree(Simplex_tree && complex_source) {
352#ifdef DEBUG_TRACES
353 std::clog << "Simplex_tree move constructor" << std::endl;
354#endif // DEBUG_TRACES
355 move_from(complex_source);
356
357 // just need to set dimension_ on source to make it available again
358 // (filtration_vect_ and members are already set from the move)
359 complex_source.dimension_ = -1;
360 }
361
364 root_members_recursive_deletion();
365 }
366
368 Simplex_tree& operator= (const Simplex_tree& complex_source) {
369#ifdef DEBUG_TRACES
370 std::clog << "Simplex_tree copy assignment" << std::endl;
371#endif // DEBUG_TRACES
372 // Self-assignment detection
373 if (&complex_source != this) {
374 // We start by deleting root_ if not empty
375 root_members_recursive_deletion();
376
377 copy_from(complex_source);
378 }
379 return *this;
380 }
381
386#ifdef DEBUG_TRACES
387 std::clog << "Simplex_tree move assignment" << std::endl;
388#endif // DEBUG_TRACES
389 // Self-assignment detection
390 if (&complex_source != this) {
391 // root_ deletion in case it was not empty
392 root_members_recursive_deletion();
393
394 move_from(complex_source);
395 }
396 return *this;
397 } // end constructor/destructor
399
400 private:
401 // Copy from complex_source to "this"
402 void copy_from(const Simplex_tree& complex_source) {
403 null_vertex_ = complex_source.null_vertex_;
404 filtration_vect_.clear();
405 dimension_ = complex_source.dimension_;
406 auto root_source = complex_source.root_;
407
408 // root members copy
409 root_.members() = Dictionary(boost::container::ordered_unique_range, root_source.members().begin(), root_source.members().end());
410 // Needs to reassign children
411 for (auto& map_el : root_.members()) {
412 map_el.second.assign_children(&root_);
413 }
414 rec_copy(&root_, &root_source);
415 }
416
418 void rec_copy(Siblings *sib, Siblings *sib_source) {
419 for (auto sh = sib->members().begin(), sh_source = sib_source->members().begin();
420 sh != sib->members().end(); ++sh, ++sh_source) {
421 if (has_children(sh_source)) {
422 Siblings * newsib = new Siblings(sib, sh_source->first);
423 newsib->members_.reserve(sh_source->second.children()->members().size());
424 for (auto & child : sh_source->second.children()->members())
425 newsib->members_.emplace_hint(newsib->members_.end(), child.first, Node(newsib, child.second.filtration()));
426 rec_copy(newsib, sh_source->second.children());
427 sh->second.assign_children(newsib);
428 }
429 }
430 }
431
432 // Move from complex_source to "this"
433 void move_from(Simplex_tree& complex_source) {
434 null_vertex_ = std::move(complex_source.null_vertex_);
435 root_ = std::move(complex_source.root_);
436 filtration_vect_ = std::move(complex_source.filtration_vect_);
437 dimension_ = complex_source.dimension_;
438
439 // Need to update root members (children->oncles and children need to point on the new root pointer)
440 for (auto& map_el : root_.members()) {
441 if (map_el.second.children() != &(complex_source.root_)) {
442 // reset children->oncles with the moved root_ pointer value
443 map_el.second.children()->oncles_ = &root_;
444 } else {
445 // if simplex is of dimension 0, oncles_ shall be nullptr
446 GUDHI_CHECK(map_el.second.children()->oncles_ == nullptr,
447 std::invalid_argument("Simplex_tree move constructor from an invalid Simplex_tree"));
448 // and children points on root_ - to be moved
449 map_el.second.assign_children(&root_);
450 }
451 }
452 }
453
454 // delete all root_.members() recursively
455 void root_members_recursive_deletion() {
456 for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) {
457 if (has_children(sh)) {
458 rec_delete(sh->second.children());
459 }
460 }
461 root_.members().clear();
462 }
463
464 // Recursive deletion
465 void rec_delete(Siblings * sib) {
466 for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
467 if (has_children(sh)) {
468 rec_delete(sh->second.children());
469 }
470 }
471 delete sib;
472 }
473
474 public:
477 if ((null_vertex_ != st2.null_vertex_) ||
478 (dimension_ != st2.dimension_ && !dimension_to_be_lowered_ && !st2.dimension_to_be_lowered_))
479 return false;
480 return rec_equal(&root_, &st2.root_);
481 }
482
485 return (!(*this == st2));
486 }
487
488 private:
490 bool rec_equal(Siblings* s1, Siblings* s2) {
491 if (s1->members().size() != s2->members().size())
492 return false;
493 for (auto sh1 = s1->members().begin(), sh2 = s2->members().begin();
494 (sh1 != s1->members().end() && sh2 != s2->members().end()); ++sh1, ++sh2) {
495 if (sh1->first != sh2->first || sh1->second.filtration() != sh2->second.filtration())
496 return false;
497 if (has_children(sh1) != has_children(sh2))
498 return false;
499 // Recursivity on children only if both have children
500 else if (has_children(sh1))
501 if (!rec_equal(sh1->second.children(), sh2->second.children()))
502 return false;
503 }
504 return true;
505 }
506
511 static Filtration_value filtration_(Simplex_handle sh) {
512 GUDHI_CHECK (sh != null_simplex(), "null simplex");
513 return sh->second.filtration();
514 }
515
516 public:
523 return sh->second.key();
524 }
525
531 return filtration_vect_[idx];
532 }
533
540 if (sh != null_simplex()) {
541 return sh->second.filtration();
542 } else {
543 return std::numeric_limits<Filtration_value>::infinity();
544 }
545 }
546
551 GUDHI_CHECK(sh != null_simplex(),
552 std::invalid_argument("Simplex_tree::assign_filtration - cannot assign filtration on null_simplex"));
553 sh->second.assign_filtration(fv);
554 }
555
561 return Dictionary_it(nullptr);
562 }
563
566 return -1;
567 }
568
572 return null_vertex_;
573 }
574
576 size_t num_vertices() const {
577 return root_.members_.size();
578 }
579
581 bool is_empty() const {
582 return root_.members_.empty();
583 }
584
585 public:
587 size_t num_simplices() {
588 return num_simplices(&root_);
589 }
590
591 private:
593 size_t num_simplices(Siblings * sib) {
594 auto sib_begin = sib->members().begin();
595 auto sib_end = sib->members().end();
596 size_t simplices_number = sib_end - sib_begin;
597 for (auto sh = sib_begin; sh != sib_end; ++sh) {
598 if (has_children(sh)) {
599 simplices_number += num_simplices(sh->second.children());
600 }
601 }
602 return simplices_number;
603 }
604
605 public:
610 Siblings * curr_sib = self_siblings(sh);
611 int dim = 0;
612 while (curr_sib != nullptr) {
613 ++dim;
614 curr_sib = curr_sib->oncles();
615 }
616 return dim - 1;
617 }
618
621 return dimension_;
622 }
623
628 int dimension() {
629 if (dimension_to_be_lowered_)
630 lower_upper_bound_dimension();
631 return dimension_;
632 }
633
636 template<class SimplexHandle>
637 bool has_children(SimplexHandle sh) const {
638 // Here we rely on the root using null_vertex(), which cannot match any real vertex.
639 return (sh->second.children()->parent() == sh->first);
640 }
641
649 template<class InputVertexRange = std::initializer_list<Vertex_handle>>
650 Simplex_handle find(const InputVertexRange & s) {
651 auto first = std::begin(s);
652 auto last = std::end(s);
653
654 if (first == last)
655 return null_simplex(); // ----->>
656
657 // Copy before sorting
658 std::vector<Vertex_handle> copy(first, last);
659 std::sort(std::begin(copy), std::end(copy));
660 return find_simplex(copy);
661 }
662
663 private:
665 Simplex_handle find_simplex(const std::vector<Vertex_handle> & simplex) {
666 Siblings * tmp_sib = &root_;
667 Dictionary_it tmp_dit;
668 auto vi = simplex.begin();
670 // Equivalent to the first iteration of the normal loop
671 GUDHI_CHECK(contiguous_vertices(), "non-contiguous vertices");
672 Vertex_handle v = *vi++;
673 if(v < 0 || v >= static_cast<Vertex_handle>(root_.members_.size()))
674 return null_simplex();
675 tmp_dit = root_.members_.begin() + v;
676 if (vi == simplex.end())
677 return tmp_dit;
678 if (!has_children(tmp_dit))
679 return null_simplex();
680 tmp_sib = tmp_dit->second.children();
681 }
682 for (;;) {
683 tmp_dit = tmp_sib->members_.find(*vi++);
684 if (tmp_dit == tmp_sib->members_.end())
685 return null_simplex();
686 if (vi == simplex.end())
687 return tmp_dit;
688 if (!has_children(tmp_dit))
689 return null_simplex();
690 tmp_sib = tmp_dit->second.children();
691 }
692 }
693
696 Simplex_handle find_vertex(Vertex_handle v) {
698 assert(contiguous_vertices());
699 return root_.members_.begin() + v;
700 } else {
701 return root_.members_.find(v);
702 }
703 }
704
705 public:
707 bool contiguous_vertices() const {
708 if (root_.members_.empty()) return true;
709 if (root_.members_.begin()->first != 0) return false;
710 if (std::prev(root_.members_.end())->first != static_cast<Vertex_handle>(root_.members_.size() - 1)) return false;
711 return true;
712 }
713
714 protected:
729 template <class RandomVertexHandleRange = std::initializer_list<Vertex_handle>>
730 std::pair<Simplex_handle, bool> insert_simplex_raw(const RandomVertexHandleRange& simplex,
732 Siblings * curr_sib = &root_;
733 std::pair<Simplex_handle, bool> res_insert;
734 auto vi = simplex.begin();
735 for (; vi != std::prev(simplex.end()); ++vi) {
736 GUDHI_CHECK(*vi != null_vertex(), "cannot use the dummy null_vertex() as a real vertex");
737 res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration));
738 if (!(has_children(res_insert.first))) {
739 res_insert.first->second.assign_children(new Siblings(curr_sib, *vi));
740 }
741 curr_sib = res_insert.first->second.children();
742 }
743 GUDHI_CHECK(*vi != null_vertex(), "cannot use the dummy null_vertex() as a real vertex");
744 res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration));
745 if (!res_insert.second) {
746 // if already in the complex
747 if (res_insert.first->second.filtration() > filtration) {
748 // if filtration value modified
749 res_insert.first->second.assign_filtration(filtration);
750 return res_insert;
751 }
752 // if filtration value unchanged
753 return std::pair<Simplex_handle, bool>(null_simplex(), false);
754 }
755 // otherwise the insertion has succeeded - size is a size_type
756 int dim = static_cast<int>(boost::size(simplex)) - 1;
757 if (dim > dimension_) {
758 // Update dimension if needed
759 dimension_ = dim;
760 }
761 return res_insert;
762 }
763
764 public:
788 template<class InputVertexRange = std::initializer_list<Vertex_handle>>
789 std::pair<Simplex_handle, bool> insert_simplex(const InputVertexRange & simplex,
791 auto first = std::begin(simplex);
792 auto last = std::end(simplex);
793
794 if (first == last)
795 return std::pair<Simplex_handle, bool>(null_simplex(), true); // ----->>
796
797 // Copy before sorting
798 std::vector<Vertex_handle> copy(first, last);
799 std::sort(std::begin(copy), std::end(copy));
800 return insert_simplex_raw(copy, filtration);
801 }
802
817 template<class InputVertexRange = std::initializer_list<Vertex_handle>>
818 std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(const InputVertexRange& Nsimplex,
820 auto first = std::begin(Nsimplex);
821 auto last = std::end(Nsimplex);
822
823 if (first == last)
824 return { null_simplex(), true }; // FIXME: false would make more sense to me.
825
826 thread_local std::vector<Vertex_handle> copy;
827 copy.clear();
828 copy.insert(copy.end(), first, last);
829 std::sort(copy.begin(), copy.end());
830 auto last_unique = std::unique(copy.begin(), copy.end());
831 copy.erase(last_unique, copy.end());
832 GUDHI_CHECK_code(
833 for (Vertex_handle v : copy)
834 GUDHI_CHECK(v != null_vertex(), "cannot use the dummy null_vertex() as a real vertex");
835 )
836 // Update dimension if needed. We could wait to see if the insertion succeeds, but I doubt there is much to gain.
837 dimension_ = (std::max)(dimension_, static_cast<int>(std::distance(copy.begin(), copy.end())) - 1);
838
839 return rec_insert_simplex_and_subfaces_sorted(root(), copy.begin(), copy.end(), filtration);
840 }
841
842 private:
843 // To insert {1,2,3,4}, we insert {2,3,4} twice, once at the root, and once below 1.
844 template<class ForwardVertexIterator>
845 std::pair<Simplex_handle, bool> rec_insert_simplex_and_subfaces_sorted(Siblings* sib,
846 ForwardVertexIterator first,
847 ForwardVertexIterator last,
848 Filtration_value filt) {
849 // An alternative strategy would be:
850 // - try to find the complete simplex, if found (and low filtration) exit
851 // - insert all the vertices at once in sib
852 // - loop over those (new or not) simplices, with a recursive call(++first, last)
853 Vertex_handle vertex_one = *first;
854 auto&& dict = sib->members();
855 auto insertion_result = dict.emplace(vertex_one, Node(sib, filt));
856 Simplex_handle simplex_one = insertion_result.first;
857 bool one_is_new = insertion_result.second;
858 if (!one_is_new) {
859 if (filtration(simplex_one) > filt) {
860 assign_filtration(simplex_one, filt);
861 } else {
862 // FIXME: this interface makes no sense, and it doesn't seem to be tested.
863 insertion_result.first = null_simplex();
864 }
865 }
866 if (++first == last) return insertion_result;
867 if (!has_children(simplex_one))
868 // TODO: have special code here, we know we are building the whole subtree from scratch.
869 simplex_one->second.assign_children(new Siblings(sib, vertex_one));
870 auto res = rec_insert_simplex_and_subfaces_sorted(simplex_one->second.children(), first, last, filt);
871 // No need to continue if the full simplex was already there with a low enough filtration value.
872 if (res.first != null_simplex()) rec_insert_simplex_and_subfaces_sorted(sib, first, last, filt);
873 return res;
874 }
875
876 public:
880 sh->second.assign_key(key);
881 }
882
886 std::pair<Simplex_handle, Simplex_handle> endpoints(Simplex_handle sh) {
887 assert(dimension(sh) == 1);
888 return { find_vertex(sh->first), find_vertex(self_siblings(sh)->parent()) };
889 }
890
892 template<class SimplexHandle>
893 static Siblings* self_siblings(SimplexHandle sh) {
894 if (sh->second.children()->parent() == sh->first)
895 return sh->second.children()->oncles();
896 else
897 return sh->second.children();
898 }
899
900 public:
903 return &root_;
904 }
905
912 void set_dimension(int dimension, bool exact=true) {
913 dimension_to_be_lowered_ = !exact;
914 dimension_ = dimension;
915 }
916
917 public:
926 filtration_vect_.clear();
927 filtration_vect_.reserve(num_simplices());
929 filtration_vect_.push_back(sh);
930
931 /* We use stable_sort here because with libstdc++ it is faster than sort.
932 * is_before_in_filtration is now a total order, but we used to call
933 * stable_sort for the following heuristic:
934 * The use of a depth-first traversal of the simplex tree, provided by
935 * complex_simplex_range(), combined with a stable sort is meant to
936 * optimize the order of simplices with same filtration value. The
937 * heuristic consists in inserting the cofaces of a simplex as soon as
938 * possible.
939 */
940#ifdef GUDHI_USE_TBB
941 tbb::parallel_sort(filtration_vect_.begin(), filtration_vect_.end(), is_before_in_filtration(this));
942#else
943 std::stable_sort(filtration_vect_.begin(), filtration_vect_.end(), is_before_in_filtration(this));
944#endif
945 }
950 if (filtration_vect_.empty()) {
952 }
953 }
959 filtration_vect_.clear();
960 }
961
962 private:
975 void rec_coface(std::vector<Vertex_handle> &vertices, Siblings *curr_sib, int curr_nbVertices,
976 std::vector<Simplex_handle>& cofaces, bool star, int nbVertices) {
977 if (!(star || curr_nbVertices <= nbVertices)) // dimension of actual simplex <= nbVertices
978 return;
979 for (Simplex_handle simplex = curr_sib->members().begin(); simplex != curr_sib->members().end(); ++simplex) {
980 if (vertices.empty()) {
981 // If we reached the end of the vertices, and the simplex has more vertices than the given simplex
982 // => we found a coface
983
984 // Add a coface if we want the star or if the number of vertices of the current simplex matches with nbVertices
985 bool addCoface = (star || curr_nbVertices == nbVertices);
986 if (addCoface)
987 cofaces.push_back(simplex);
988 if ((!addCoface || star) && has_children(simplex)) // Rec call
989 rec_coface(vertices, simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
990 } else {
991 if (simplex->first == vertices.back()) {
992 // If curr_sib matches with the top vertex
993 bool equalDim = (star || curr_nbVertices == nbVertices); // dimension of actual simplex == nbVertices
994 bool addCoface = vertices.size() == 1 && equalDim;
995 if (addCoface)
996 cofaces.push_back(simplex);
997 if ((!addCoface || star) && has_children(simplex)) {
998 // Rec call
999 Vertex_handle tmp = vertices.back();
1000 vertices.pop_back();
1001 rec_coface(vertices, simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
1002 vertices.push_back(tmp);
1003 }
1004 } else if (simplex->first > vertices.back()) {
1005 return;
1006 } else {
1007 // (simplex->first < vertices.back()
1008 if (has_children(simplex))
1009 rec_coface(vertices, simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
1010 }
1011 }
1012 }
1013 }
1014
1015 public:
1022 return cofaces_simplex_range(simplex, 0);
1023 }
1024
1033 Cofaces_simplex_range cofaces;
1034 // codimension must be positive or null integer
1035 assert(codimension >= 0);
1037 std::vector<Vertex_handle> copy(rg.begin(), rg.end());
1038 if (codimension + static_cast<int>(copy.size()) > dimension_ + 1 ||
1039 (codimension == 0 && static_cast<int>(copy.size()) > dimension_)) // n+codimension greater than dimension_
1040 return cofaces;
1041 // must be sorted in decreasing order
1042 assert(std::is_sorted(copy.begin(), copy.end(), std::greater<Vertex_handle>()));
1043 bool star = codimension == 0;
1044 rec_coface(copy, &root_, 1, cofaces, star, codimension + static_cast<int>(copy.size()));
1045 return cofaces;
1046 }
1047
1048 private:
1056 bool reverse_lexicographic_order(Simplex_handle sh1, Simplex_handle sh2) {
1059 Simplex_vertex_iterator it1 = rg1.begin();
1060 Simplex_vertex_iterator it2 = rg2.begin();
1061 while (it1 != rg1.end() && it2 != rg2.end()) {
1062 if (*it1 == *it2) {
1063 ++it1;
1064 ++it2;
1065 } else {
1066 return *it1 < *it2;
1067 }
1068 }
1069 return ((it1 == rg1.end()) && (it2 != rg2.end()));
1070 }
1071
1078 struct is_before_in_filtration {
1079 explicit is_before_in_filtration(Simplex_tree * st)
1080 : st_(st) { }
1081
1082 bool operator()(const Simplex_handle sh1, const Simplex_handle sh2) const {
1083 // Not using st_->filtration(sh1) because it uselessly tests for null_simplex.
1084 if (sh1->second.filtration() != sh2->second.filtration()) {
1085 return sh1->second.filtration() < sh2->second.filtration();
1086 }
1087 // is sh1 a proper subface of sh2
1088 return st_->reverse_lexicographic_order(sh1, sh2);
1089 }
1090
1091 Simplex_tree * st_;
1092 };
1093
1094 public:
1118 template<class OneSkeletonGraph>
1119 void insert_graph(const OneSkeletonGraph& skel_graph) {
1120 // the simplex tree must be empty
1121 assert(num_simplices() == 0);
1122
1123 // is there a better way to let the compiler know that we don't mean Simplex_tree::num_vertices?
1124 using boost::num_vertices;
1125
1126 if (num_vertices(skel_graph) == 0) {
1127 return;
1128 }
1129 if (num_edges(skel_graph) == 0) {
1130 dimension_ = 0;
1131 } else {
1132 dimension_ = 1;
1133 }
1134
1135 root_.members_.reserve(num_vertices(skel_graph)); // probably useless in most cases
1136 auto verts = vertices(skel_graph) | boost::adaptors::transformed([&](auto v){
1137 return Dit_value_t(v, Node(&root_, get(vertex_filtration_t(), skel_graph, v))); });
1138 root_.members_.insert(boost::begin(verts), boost::end(verts));
1139 // This automatically sorts the vertices, the graph concept doesn't guarantee the order in which we iterate.
1140
1141 std::pair<typename boost::graph_traits<OneSkeletonGraph>::edge_iterator,
1142 typename boost::graph_traits<OneSkeletonGraph>::edge_iterator> boost_edges = edges(skel_graph);
1143 // boost_edges.first is the equivalent to boost_edges.begin()
1144 // boost_edges.second is the equivalent to boost_edges.end()
1145 for (; boost_edges.first != boost_edges.second; boost_edges.first++) {
1146 auto edge = *(boost_edges.first);
1147 auto u = source(edge, skel_graph);
1148 auto v = target(edge, skel_graph);
1149 if (u == v) throw std::invalid_argument("Self-loops are not simplicial");
1150 // We cannot skip edges with the wrong orientation and expect them to
1151 // come a second time with the right orientation, that does not always
1152 // happen in practice. emplace() should be a NOP when an element with the
1153 // same key is already there, so seeing the same edge multiple times is
1154 // ok.
1155 // Should we actually forbid multiple edges? That would be consistent
1156 // with rejecting self-loops.
1157 if (v < u) std::swap(u, v);
1158 auto sh = find_vertex(u);
1159 if (!has_children(sh)) {
1160 sh->second.assign_children(new Siblings(&root_, sh->first));
1161 }
1162
1163 sh->second.children()->members().emplace(v,
1164 Node(sh->second.children(), get(edge_filtration_t(), skel_graph, edge)));
1165 }
1166 }
1167
1175 template <class VertexRange>
1176 void insert_batch_vertices(VertexRange const& vertices, Filtration_value filt = 0) {
1177 auto verts = vertices | boost::adaptors::transformed([&](auto v){
1178 return Dit_value_t(v, Node(&root_, filt)); });
1179 root_.members_.insert(boost::begin(verts), boost::end(verts));
1180 if (dimension_ < 0 && !root_.members_.empty()) dimension_ = 0;
1181 }
1182
1194 void expansion(int max_dim) {
1195 if (max_dim <= 1) return;
1196 clear_filtration(); // Drop the cache.
1197 dimension_ = max_dim;
1198 for (Dictionary_it root_it = root_.members_.begin();
1199 root_it != root_.members_.end(); ++root_it) {
1200 if (has_children(root_it)) {
1201 siblings_expansion(root_it->second.children(), max_dim - 1);
1202 }
1203 }
1204 dimension_ = max_dim - dimension_;
1205 }
1206
1207 private:
1209 void siblings_expansion(Siblings * siblings, // must contain elements
1210 int k) {
1211 if (dimension_ > k) {
1212 dimension_ = k;
1213 }
1214 if (k == 0)
1215 return;
1216 Dictionary_it next = siblings->members().begin();
1217 ++next;
1218
1219 thread_local std::vector<std::pair<Vertex_handle, Node> > inter;
1220 for (Dictionary_it s_h = siblings->members().begin();
1221 s_h != siblings->members().end(); ++s_h, ++next) {
1222 Simplex_handle root_sh = find_vertex(s_h->first);
1223 if (has_children(root_sh)) {
1224 intersection(
1225 inter, // output intersection
1226 next, // begin
1227 siblings->members().end(), // end
1228 root_sh->second.children()->members().begin(),
1229 root_sh->second.children()->members().end(),
1230 s_h->second.filtration());
1231 if (inter.size() != 0) {
1232 Siblings * new_sib = new Siblings(siblings, // oncles
1233 s_h->first, // parent
1234 inter); // boost::container::ordered_unique_range_t
1235 inter.clear();
1236 s_h->second.assign_children(new_sib);
1237 siblings_expansion(new_sib, k - 1);
1238 } else {
1239 // ensure the children property
1240 s_h->second.assign_children(siblings);
1241 inter.clear();
1242 }
1243 }
1244 }
1245 }
1246
1249 static void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection,
1250 Dictionary_it begin1, Dictionary_it end1,
1251 Dictionary_it begin2, Dictionary_it end2,
1252 Filtration_value filtration_) {
1253 if (begin1 == end1 || begin2 == end2)
1254 return; // ----->>
1255 while (true) {
1256 if (begin1->first == begin2->first) {
1257 Filtration_value filt = (std::max)({begin1->second.filtration(), begin2->second.filtration(), filtration_});
1258 intersection.emplace_back(begin1->first, Node(nullptr, filt));
1259 if (++begin1 == end1 || ++begin2 == end2)
1260 return; // ----->>
1261 } else if (begin1->first < begin2->first) {
1262 if (++begin1 == end1)
1263 return;
1264 } else /* begin1->first > begin2->first */ {
1265 if (++begin2 == end2)
1266 return; // ----->>
1267 }
1268 }
1269 }
1270
1271 public:
1290 template< typename Blocker >
1291 void expansion_with_blockers(int max_dim, Blocker block_simplex) {
1292 // Loop must be from the end to the beginning, as higher dimension simplex are always on the left part of the tree
1293 for (auto& simplex : boost::adaptors::reverse(root_.members())) {
1294 if (has_children(&simplex)) {
1295 siblings_expansion_with_blockers(simplex.second.children(), max_dim, max_dim - 1, block_simplex);
1296 }
1297 }
1298 }
1299
1300 private:
1302 template< typename Blocker >
1303 void siblings_expansion_with_blockers(Siblings* siblings, int max_dim, int k, Blocker block_simplex) {
1304 if (dimension_ < max_dim - k) {
1305 dimension_ = max_dim - k;
1306 }
1307 if (k == 0)
1308 return;
1309 // No need to go deeper
1310 if (siblings->members().size() < 2)
1311 return;
1312 // Reverse loop starting before the last one for 'next' to be the last one
1313 for (auto simplex = siblings->members().rbegin() + 1; simplex != siblings->members().rend(); simplex++) {
1314 std::vector<std::pair<Vertex_handle, Node> > intersection;
1315 for(auto next = siblings->members().rbegin(); next != simplex; next++) {
1316 bool to_be_inserted = true;
1317 Filtration_value filt = simplex->second.filtration();
1318 // If all the boundaries are present, 'next' needs to be inserted
1319 for (Simplex_handle border : boundary_simplex_range(simplex)) {
1320 Simplex_handle border_child = find_child(border, next->first);
1321 if (border_child == null_simplex()) {
1322 to_be_inserted=false;
1323 break;
1324 }
1325 filt = (std::max)(filt, filtration(border_child));
1326 }
1327 if (to_be_inserted) {
1328 intersection.emplace_back(next->first, Node(nullptr, filt));
1329 }
1330 }
1331 if (intersection.size() != 0) {
1332 // Reverse the order to insert
1333 Siblings * new_sib = new Siblings(siblings, // oncles
1334 simplex->first, // parent
1335 boost::adaptors::reverse(intersection)); // boost::container::ordered_unique_range_t
1336 simplex->second.assign_children(new_sib);
1337 std::vector<Vertex_handle> blocked_new_sib_vertex_list;
1338 // As all intersections are inserted, we can call the blocker function on all new_sib members
1339 for (auto new_sib_member = new_sib->members().begin();
1340 new_sib_member != new_sib->members().end();
1341 new_sib_member++) {
1342 bool blocker_result = block_simplex(new_sib_member);
1343 // new_sib member has been blocked by the blocker function
1344 // add it to the list to be removed - do not perform it while looping on it
1345 if (blocker_result) {
1346 blocked_new_sib_vertex_list.push_back(new_sib_member->first);
1347 }
1348 }
1349 if (blocked_new_sib_vertex_list.size() == new_sib->members().size()) {
1350 // Specific case where all have to be deleted
1351 delete new_sib;
1352 // ensure the children property
1353 simplex->second.assign_children(siblings);
1354 } else {
1355 for (auto& blocked_new_sib_member : blocked_new_sib_vertex_list) {
1356 new_sib->members().erase(blocked_new_sib_member);
1357 }
1358 // ensure recursive call
1359 siblings_expansion_with_blockers(new_sib, max_dim, k - 1, block_simplex);
1360 }
1361 } else {
1362 // ensure the children property
1363 simplex->second.assign_children(siblings);
1364 }
1365 }
1366 }
1367
1372 Simplex_handle find_child(Simplex_handle sh, Vertex_handle vh) const {
1373 if (!has_children(sh))
1374 return null_simplex();
1375
1376 Simplex_handle child = sh->second.children()->find(vh);
1377 // Specific case of boost::flat_map does not find, returns boost::flat_map::end()
1378 // in simplex tree we want a null_simplex()
1379 if (child == sh->second.children()->members().end())
1380 return null_simplex();
1381
1382 return child;
1383 }
1384
1385 public:
1392 void print_hasse(std::ostream& os) {
1393 os << num_simplices() << " " << std::endl;
1394 for (auto sh : filtration_simplex_range()) {
1395 os << dimension(sh) << " ";
1396 for (auto b_sh : boundary_simplex_range(sh)) {
1397 os << key(b_sh) << " ";
1398 }
1399 os << filtration(sh) << " \n";
1400 }
1401 }
1402
1403 public:
1411 bool modified = false;
1412 // Loop must be from the end to the beginning, as higher dimension simplex are always on the left part of the tree
1413 for (auto& simplex : boost::adaptors::reverse(root_.members())) {
1414 if (has_children(&simplex)) {
1415 modified |= rec_make_filtration_non_decreasing(simplex.second.children());
1416 }
1417 }
1418 if(modified)
1419 clear_filtration(); // Drop the cache.
1420 return modified;
1421 }
1422
1423 private:
1428 bool rec_make_filtration_non_decreasing(Siblings * sib) {
1429 bool modified = false;
1430
1431 // Loop must be from the end to the beginning, as higher dimension simplex are always on the left part of the tree
1432 for (auto& simplex : boost::adaptors::reverse(sib->members())) {
1433 // Find the maximum filtration value in the border
1435 Boundary_simplex_iterator max_border = std::max_element(std::begin(boundary), std::end(boundary),
1436 [](Simplex_handle sh1, Simplex_handle sh2) {
1437 return filtration(sh1) < filtration(sh2);
1438 });
1439
1440 Filtration_value max_filt_border_value = filtration(*max_border);
1441 // Replacing if(f<max) with if(!(f>=max)) would mean that if f is NaN, we replace it with the max of the children.
1442 // That seems more useful than keeping NaN.
1443 if (!(simplex.second.filtration() >= max_filt_border_value)) {
1444 // Store the filtration modification information
1445 modified = true;
1446 simplex.second.assign_filtration(max_filt_border_value);
1447 }
1448 if (has_children(&simplex)) {
1449 modified |= rec_make_filtration_non_decreasing(simplex.second.children());
1450 }
1451 }
1452 // Make the modified information to be traced by upper call
1453 return modified;
1454 }
1455
1456 public:
1465 bool modified = rec_prune_above_filtration(root(), filtration);
1466 if(modified)
1467 clear_filtration(); // Drop the cache.
1468 return modified;
1469 }
1470
1471 private:
1472 bool rec_prune_above_filtration(Siblings* sib, Filtration_value filt) {
1473 auto&& list = sib->members();
1474 auto last = std::remove_if(list.begin(), list.end(), [this,filt](Dit_value_t& simplex) {
1475 if (simplex.second.filtration() <= filt) return false;
1476 if (has_children(&simplex)) rec_delete(simplex.second.children());
1477 // dimension may need to be lowered
1478 dimension_to_be_lowered_ = true;
1479 return true;
1480 });
1481
1482 bool modified = (last != list.end());
1483 if (last == list.begin() && sib != root()) {
1484 // Removing the whole siblings, parent becomes a leaf.
1485 sib->oncles()->members()[sib->parent()].assign_children(sib->oncles());
1486 delete sib;
1487 // dimension may need to be lowered
1488 dimension_to_be_lowered_ = true;
1489 return true;
1490 } else {
1491 // Keeping some elements of siblings. Remove the others, and recurse in the remaining ones.
1492 list.erase(last, list.end());
1493 for (auto&& simplex : list)
1494 if (has_children(&simplex))
1495 modified |= rec_prune_above_filtration(simplex.second.children(), filt);
1496 }
1497 return modified;
1498 }
1499
1500 public:
1506 if (dimension >= dimension_)
1507 return false;
1508
1509 bool modified = false;
1510 if (dimension < 0) {
1511 if (num_vertices() > 0) {
1512 root_members_recursive_deletion();
1513 modified = true;
1514 }
1515 // Force dimension to -1, in case user calls `prune_above_dimension(-10)`
1516 dimension = -1;
1517 } else {
1518 modified = rec_prune_above_dimension(root(), dimension, 0);
1519 }
1520 if(modified) {
1521 // Thanks to `if (dimension >= dimension_)` and dimension forced to -1 `if (dimension < 0)`, we know the new dimension
1522 dimension_ = dimension;
1523 clear_filtration(); // Drop the cache.
1524 }
1525 return modified;
1526 }
1527
1528 private:
1529 bool rec_prune_above_dimension(Siblings* sib, int dim, int actual_dim) {
1530 bool modified = false;
1531 auto&& list = sib->members();
1532
1533 for (auto&& simplex : list)
1534 if (has_children(&simplex)) {
1535 if (actual_dim >= dim) {
1536 rec_delete(simplex.second.children());
1537 simplex.second.assign_children(sib);
1538 modified = true;
1539 } else {
1540 modified |= rec_prune_above_dimension(simplex.second.children(), dim, actual_dim + 1);
1541 }
1542 }
1543
1544 return modified;
1545 }
1546
1547 private:
1553 bool lower_upper_bound_dimension() {
1554 // reset automatic detection to recompute
1555 dimension_to_be_lowered_ = false;
1556 int new_dimension = -1;
1557 // Browse the tree from the left to the right as higher dimension cells are more likely on the left part of the tree
1558 for (Simplex_handle sh : complex_simplex_range()) {
1559#ifdef DEBUG_TRACES
1560 for (auto vertex : simplex_vertex_range(sh)) {
1561 std::clog << " " << vertex;
1562 }
1563 std::clog << std::endl;
1564#endif // DEBUG_TRACES
1565
1566 int sh_dimension = dimension(sh);
1567 if (sh_dimension >= dimension_)
1568 // Stop browsing as soon as the dimension is reached, no need to go further
1569 return false;
1570 new_dimension = (std::max)(new_dimension, sh_dimension);
1571 }
1572 dimension_ = new_dimension;
1573 return true;
1574 }
1575
1576
1577 public:
1587 // Guarantee the simplex has no children
1588 GUDHI_CHECK(!has_children(sh),
1589 std::invalid_argument("Simplex_tree::remove_maximal_simplex - argument has children"));
1590
1591 // Simplex is a leaf, it means the child is the Siblings owning the leaf
1592 Siblings* child = sh->second.children();
1593
1594 if ((child->size() > 1) || (child == root())) {
1595 // Not alone, just remove it from members
1596 // Special case when child is the root of the simplex tree, just remove it from members
1597 child->erase(sh);
1598 } else {
1599 // Sibling is emptied : must be deleted, and its parent must point on his own Sibling
1600 child->oncles()->members().at(child->parent()).assign_children(child->oncles());
1601 delete child;
1602 // dimension may need to be lowered
1603 dimension_to_be_lowered_ = true;
1604 }
1605 }
1606
1623 std::pair<Filtration_value, Extended_simplex_type> decode_extended_filtration(Filtration_value f, const Extended_filtration_data& efd){
1624 std::pair<Filtration_value, Extended_simplex_type> p;
1625 Filtration_value minval = efd.minval;
1626 Filtration_value maxval = efd.maxval;
1627 if (f >= -2 && f <= -1){
1628 p.first = minval + (maxval-minval)*(f + 2); p.second = Extended_simplex_type::UP;
1629 }
1630 else if (f >= 1 && f <= 2){
1631 p.first = minval - (maxval-minval)*(f - 2); p.second = Extended_simplex_type::DOWN;
1632 }
1633 else{
1634 p.first = std::numeric_limits<Filtration_value>::quiet_NaN(); p.second = Extended_simplex_type::EXTRA;
1635 }
1636 return p;
1637 };
1638
1652 Extended_filtration_data extend_filtration() {
1653 clear_filtration(); // Drop the cache.
1654
1655 // Compute maximum and minimum of filtration values
1656 Vertex_handle maxvert = std::numeric_limits<Vertex_handle>::min();
1657 Filtration_value minval = std::numeric_limits<Filtration_value>::infinity();
1658 Filtration_value maxval = -std::numeric_limits<Filtration_value>::infinity();
1659 for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh){
1660 Filtration_value f = this->filtration(sh);
1661 minval = std::min(minval, f);
1662 maxval = std::max(maxval, f);
1663 maxvert = std::max(sh->first, maxvert);
1664 }
1665
1666 GUDHI_CHECK(maxvert < std::numeric_limits<Vertex_handle>::max(), std::invalid_argument("Simplex_tree contains a vertex with the largest Vertex_handle"));
1667 maxvert += 1;
1668
1669 Simplex_tree st_copy = *this;
1670
1671 // Add point for coning the simplicial complex
1672 this->insert_simplex_raw({maxvert}, -3);
1673
1674 // For each simplex
1675 std::vector<Vertex_handle> vr;
1676 for (auto sh_copy : st_copy.complex_simplex_range()){
1677
1678 // Locate simplex
1679 vr.clear();
1680 for (auto vh : st_copy.simplex_vertex_range(sh_copy)){
1681 vr.push_back(vh);
1682 }
1683 auto sh = this->find(vr);
1684
1685 // Create cone on simplex
1686 vr.push_back(maxvert);
1687 if (this->dimension(sh) == 0){
1688 Filtration_value v = this->filtration(sh);
1689 Filtration_value scaled_v = (v-minval)/(maxval-minval);
1690 // Assign ascending value between -2 and -1 to vertex
1691 this->assign_filtration(sh, -2 + scaled_v);
1692 // Assign descending value between 1 and 2 to cone on vertex
1693 this->insert_simplex(vr, 2 - scaled_v);
1694 }
1695 else{
1696 // Assign value -3 to simplex and cone on simplex
1697 this->assign_filtration(sh, -3);
1698 this->insert_simplex(vr, -3);
1699 }
1700 }
1701
1702 // Automatically assign good values for simplices
1704
1705 // Return the filtration data
1706 Extended_filtration_data efd(minval, maxval);
1707 return efd;
1708 }
1709
1715 auto filt = filtration_(sh);
1716 for(auto v : simplex_vertex_range(sh))
1717 if(filtration_(find_vertex(v)) == filt)
1718 return v;
1719 return null_vertex();
1720 }
1721
1729 // See issue #251 for potential speed improvements.
1730 auto&& vertices = simplex_vertex_range(sh); // vertices in decreasing order
1731 auto end = std::end(vertices);
1732 auto vi = std::begin(vertices);
1733 GUDHI_CHECK(vi != end, "empty simplex");
1734 auto v0 = *vi;
1735 ++vi;
1736 GUDHI_CHECK(vi != end, "simplex of dimension 0");
1737 if(std::next(vi) == end) return sh; // shortcut for dimension 1
1738 boost::container::static_vector<Vertex_handle, 40> suffix;
1739 suffix.push_back(v0);
1740 auto filt = filtration_(sh);
1741 do
1742 {
1743 Vertex_handle v = *vi;
1744 auto&& children1 = find_vertex(v)->second.children()->members_;
1745 for(auto w : suffix){
1746 // Can we take advantage of the fact that suffix is ordered?
1747 Simplex_handle s = children1.find(w);
1748 if(filtration_(s) == filt)
1749 return s;
1750 }
1751 suffix.push_back(v);
1752 }
1753 while(++vi != end);
1754 return null_simplex();
1755 }
1756
1762 auto filt = filtration_(sh);
1763 // Naive implementation, it can be sped up.
1764 for(auto b : boundary_simplex_range(sh))
1765 if(filtration_(b) == filt)
1767 return sh; // None of its faces has the same filtration.
1768 }
1769
1770 public:
1779 void reset_filtration(Filtration_value filt_value, int min_dim = 0) {
1780 rec_reset_filtration(&root_, filt_value, min_dim);
1781 clear_filtration(); // Drop the cache.
1782 }
1783
1784 private:
1790 void rec_reset_filtration(Siblings * sib, Filtration_value filt_value, int min_depth) {
1791 for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
1792 if (min_depth <= 0) {
1793 sh->second.assign_filtration(filt_value);
1794 }
1795 if (has_children(sh)) {
1796 rec_reset_filtration(sh->second.children(), filt_value, min_depth - 1);
1797 }
1798 }
1799 }
1800
1801 public:
1809 std::size_t get_serialization_size() {
1810 const std::size_t vh_byte_size = sizeof(Vertex_handle);
1811 const std::size_t fv_byte_size = SimplexTreeOptions::store_filtration ? sizeof(Filtration_value) : 0;
1812 const std::size_t buffer_byte_size = vh_byte_size + num_simplices() * (fv_byte_size + 2 * vh_byte_size);
1813#ifdef DEBUG_TRACES
1814 std::clog << "Gudhi::simplex_tree::get_serialization_size - buffer size = " << buffer_byte_size << std::endl;
1815#endif // DEBUG_TRACES
1816 return buffer_byte_size;
1817 }
1818
1829 /* Let's take the following simplicial complex as example: */
1830 /* (vertices are represented as letters to ease the understanding) */
1831 /* o---o---o */
1832 /* a b\X/c */
1833 /* o */
1834 /* d */
1835 /* The simplex tree is: */
1836 /* a o b o c o d o */
1837 /* | |\ | */
1838 /* b o c o o d o d */
1839 /* | */
1840 /* d o */
1841 /* The serialization is (without filtration values that comes right after vertex handle value): */
1842 /* 04(number of vertices)0a 0b 0c 0d(list of vertices)01(number of [a] children)0b([a,b] simplex) */
1843 /* 00(number of [a,b] children)02(number of [b] children)0c 0d(list of [b] children)01(number of [b,c] children) */
1844 /* 0d(list of [b,c] children)00(number of [b,c,d] children)00(number of [b,d] children)01(number of [c] children) */
1845 /* 0d(list of [c] children)00(number of [b,d] children)00(number of [d] children) */
1846 /* Without explanation and with filtration values: */
1847 /* 04 0a F(a) 0b F(b) 0c F(c) 0d F(d) 01 0b F(a,b) 00 02 0c F(b,c) 0d F(b,d) 01 0d F(b,c,d) 00 00 01 0d F(c,d) 00 00 */
1848 void serialize(char* buffer, const std::size_t buffer_size) {
1849 char* buffer_end = rec_serialize(&root_, buffer);
1850 if (static_cast<std::size_t>(buffer_end - buffer) != buffer_size)
1851 throw std::invalid_argument("Serialization does not match end of buffer");
1852 }
1853
1854 private:
1856 char* rec_serialize(Siblings *sib, char* buffer) {
1857 char* ptr = buffer;
1858 ptr = Gudhi::simplex_tree::serialize_trivial(static_cast<Vertex_handle>(sib->members().size()), ptr);
1859#ifdef DEBUG_TRACES
1860 std::clog << "\n" << sib->members().size() << " : ";
1861#endif // DEBUG_TRACES
1862 for (auto& map_el : sib->members()) {
1863 ptr = Gudhi::simplex_tree::serialize_trivial(map_el.first, ptr); // Vertex
1865 ptr = Gudhi::simplex_tree::serialize_trivial(map_el.second.filtration(), ptr); // Filtration
1866#ifdef DEBUG_TRACES
1867 std::clog << " [ " << map_el.first << " | " << map_el.second.filtration() << " ] ";
1868#endif // DEBUG_TRACES
1869 }
1870 for (auto& map_el : sib->members()) {
1871 if (has_children(&map_el)) {
1872 ptr = rec_serialize(map_el.second.children(), ptr);
1873 } else {
1874 ptr = Gudhi::simplex_tree::serialize_trivial(static_cast<Vertex_handle>(0), ptr);
1875#ifdef DEBUG_TRACES
1876 std::cout << "\n0 : ";
1877#endif // DEBUG_TRACES
1878 }
1879 }
1880 return ptr;
1881 }
1882
1883 public:
1897 void deserialize(const char* buffer, const std::size_t buffer_size) {
1898 GUDHI_CHECK(num_vertices() == 0, std::logic_error("Simplex_tree::deserialize - Simplex_tree must be empty"));
1899 const char* ptr = buffer;
1900 // Needs to read size before recursivity to manage new siblings for children
1901 Vertex_handle members_size;
1902 ptr = Gudhi::simplex_tree::deserialize_trivial(members_size, ptr);
1903 ptr = rec_deserialize(&root_, members_size, ptr, 0);
1904 if (static_cast<std::size_t>(ptr - buffer) != buffer_size) {
1905 throw std::invalid_argument("Deserialization does not match end of buffer");
1906 }
1907 }
1908
1909 private:
1911 const char* rec_deserialize(Siblings *sib, Vertex_handle members_size, const char* ptr, int dim) {
1912 // In case buffer is just a 0 char
1913 if (members_size > 0) {
1914 sib->members_.reserve(members_size);
1915 Vertex_handle vertex;
1917 for (Vertex_handle idx = 0; idx < members_size; idx++) {
1918 ptr = Gudhi::simplex_tree::deserialize_trivial(vertex, ptr);
1920 ptr = Gudhi::simplex_tree::deserialize_trivial(filtration, ptr);
1921 // Default is no children
1922 sib->members_.emplace_hint(sib->members_.end(), vertex, Node(sib, filtration));
1923 } else {
1924 // Default is no children
1925 sib->members_.emplace_hint(sib->members_.end(), vertex, Node(sib));
1926 }
1927 }
1928 Vertex_handle child_size;
1929 for (auto& map_el : sib->members()) {
1930 ptr = Gudhi::simplex_tree::deserialize_trivial(child_size, ptr);
1931 if (child_size > 0) {
1932 Siblings* child = new Siblings(sib, map_el.first);
1933 map_el.second.assign_children(child);
1934 ptr = rec_deserialize(child, child_size, ptr, dim + 1);
1935 }
1936 }
1937 if (dim > dimension_) {
1938 // Update dimension if needed
1939 dimension_ = dim;
1940 }
1941 }
1942 return ptr;
1943 }
1944
1945 private:
1946 Vertex_handle null_vertex_;
1949 Siblings root_;
1951 std::vector<Simplex_handle> filtration_vect_;
1953 int dimension_;
1954 bool dimension_to_be_lowered_ = false;
1955};
1956
1957// Print a Simplex_tree in os.
1958template<typename...T>
1959std::ostream& operator<<(std::ostream & os, Simplex_tree<T...> & st) {
1960 for (auto sh : st.filtration_simplex_range()) {
1961 os << st.dimension(sh) << " ";
1962 for (auto v : st.simplex_vertex_range(sh)) {
1963 os << v << " ";
1964 }
1965 os << st.filtration(sh) << "\n"; // TODO(VR): why adding the key ?? not read ?? << " " << st.key(sh) << " \n";
1966 }
1967 return os;
1968}
1969
1970template<typename...T>
1971std::istream& operator>>(std::istream & is, Simplex_tree<T...> & st) {
1972 typedef Simplex_tree<T...> ST;
1973 std::vector<typename ST::Vertex_handle> simplex;
1974 typename ST::Filtration_value fil;
1975 int max_dim = -1;
1976 while (read_simplex(is, simplex, fil)) {
1977 // read all simplices in the file as a list of vertices
1978 // Warning : simplex_size needs to be casted in int - Can be 0
1979 int dim = static_cast<int> (simplex.size() - 1);
1980 if (max_dim < dim) {
1981 max_dim = dim;
1982 }
1983 // insert every simplex in the simplex tree
1984 st.insert_simplex(simplex, fil);
1985 simplex.clear();
1986 }
1987 st.set_dimension(max_dim);
1988
1989 return is;
1990}
1991
1998 typedef int Vertex_handle;
1999 typedef double Filtration_value;
2000 typedef std::uint32_t Simplex_key;
2001 static const bool store_key = true;
2002 static const bool store_filtration = true;
2003 static const bool contiguous_vertices = false;
2004};
2005
2014 typedef int Vertex_handle;
2015 typedef float Filtration_value;
2016 typedef std::uint32_t Simplex_key;
2017 static const bool store_key = true;
2018 static const bool store_filtration = true;
2019 static const bool contiguous_vertices = true;
2020};
2021 // end addtogroup simplex_tree
2023
2024} // namespace Gudhi
2025
2026#endif // SIMPLEX_TREE_H_
Extended simplex type data structure for representing the type of simplices in an extended filtration...
Iterator over the simplices of the boundary of a simplex and their opposite vertices.
Definition: Simplex_tree_iterators.h:190
Iterator over the simplices of the boundary of a simplex.
Definition: Simplex_tree_iterators.h:83
Iterator over the simplices of a simplicial complex.
Definition: Simplex_tree_iterators.h:300
Data structure to store a set of nodes in a SimplexTree sharing the same parent node.
Definition: Simplex_tree_siblings.h:31
Iterator over the vertices of a simplex in a SimplexTree.
Definition: Simplex_tree_iterators.h:38
Iterator over the simplices of the skeleton of a given dimension of the simplicial complex.
Definition: Simplex_tree_iterators.h:374
Simplex Tree data structure for representing simplicial complexes.
Definition: Simplex_tree.h:83
static Siblings * self_siblings(SimplexHandle sh)
Definition: Simplex_tree.h:893
Simplex_tree(Simplex_tree &&complex_source)
User-defined move constructor relocates the whole tree structure.
Definition: Simplex_tree.h:351
bool operator==(Simplex_tree &st2)
Checks if two simplex trees are equal.
Definition: Simplex_tree.h:476
Options::Filtration_value Filtration_value
Type for the value of the filtration function.
Definition: Simplex_tree.h:90
Cofaces_simplex_range cofaces_simplex_range(const Simplex_handle simplex, int codimension)
Compute the cofaces of a n simplex.
Definition: Simplex_tree.h:1032
Simplex_tree_boundary_opposite_vertex_simplex_iterator< Simplex_tree > Boundary_opposite_vertex_simplex_iterator
Iterator over the simplices of the boundary of a simplex and their opposite vertices.
Definition: Simplex_tree.h:200
void reset_filtration(Filtration_value filt_value, int min_dim=0)
This function resets the filtration value of all the simplices of dimension at least min_dim....
Definition: Simplex_tree.h:1779
void assign_filtration(Simplex_handle sh, Filtration_value fv)
Sets the filtration value of a simplex.
Definition: Simplex_tree.h:550
bool make_filtration_non_decreasing()
This function ensures that each simplex has a higher filtration value than its faces by increasing th...
Definition: Simplex_tree.h:1410
Dictionary::iterator Simplex_handle
Handle type to a simplex contained in the simplicial complex represented by the simplex tree.
Definition: Simplex_tree.h:156
Vertex_handle vertex_with_same_filtration(Simplex_handle sh)
Returns a vertex of sh that has the same filtration value as sh if it exists, and null_vertex() other...
Definition: Simplex_tree.h:1714
Filtration_simplex_range const & filtration_simplex_range(Indexing_tag=Indexing_tag())
Returns a range over the simplices of the simplicial complex, in the order of the filtration.
Definition: Simplex_tree.h:275
bool prune_above_dimension(int dimension)
Remove all simplices of dimension greater than a given value.
Definition: Simplex_tree.h:1505
std::pair< Simplex_handle, bool > insert_simplex(const InputVertexRange &simplex, Filtration_value filtration=0)
Insert a simplex, represented by a range of Vertex_handles, in the simplicial complex.
Definition: Simplex_tree.h:789
Vertex_handle null_vertex() const
Returns a Vertex_handle different from all Vertex_handles associated to the vertices of the simplicia...
Definition: Simplex_tree.h:571
boost::iterator_range< Boundary_simplex_iterator > Boundary_simplex_range
Range over the simplices of the boundary of a simplex.
Definition: Simplex_tree.h:196
Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) const
Returns a range over the vertices of a simplex.
Definition: Simplex_tree.h:286
void clear_filtration()
Clears the filtration cache produced by initialize_filtration().
Definition: Simplex_tree.h:958
Simplex_tree_simplex_vertex_iterator< Simplex_tree > Simplex_vertex_iterator
Iterator over the vertices of a simplex.
Definition: Simplex_tree.h:186
static Simplex_key key(Simplex_handle sh)
Returns the key associated to a simplex.
Definition: Simplex_tree.h:522
static Filtration_value filtration(Simplex_handle sh)
Returns the filtration value of a simplex.
Definition: Simplex_tree.h:539
bool operator!=(Simplex_tree &st2)
Checks if two simplex trees are different.
Definition: Simplex_tree.h:484
bool has_children(SimplexHandle sh) const
Returns true if the node in the simplex tree pointed by sh has children.
Definition: Simplex_tree.h:637
Cofaces_simplex_range star_simplex_range(const Simplex_handle simplex)
Compute the star of a n simplex.
Definition: Simplex_tree.h:1021
void expansion_with_blockers(int max_dim, Blocker block_simplex)
Expands a simplex tree containing only a graph. Simplices corresponding to cliques in the graph are a...
Definition: Simplex_tree.h:1291
boost::iterator_range< Simplex_vertex_iterator > Simplex_vertex_range
Range over the vertices of a simplex.
Definition: Simplex_tree.h:188
void remove_maximal_simplex(Simplex_handle sh)
Remove a maximal simplex.
Definition: Simplex_tree.h:1586
std::pair< Simplex_handle, Simplex_handle > endpoints(Simplex_handle sh)
Definition: Simplex_tree.h:886
void assign_key(Simplex_handle sh, Simplex_key key)
Assign a value 'key' to the key of the simplex represented by the Simplex_handle 'sh'.
Definition: Simplex_tree.h:879
Options::Simplex_key Simplex_key
Key associated to each simplex.
Definition: Simplex_tree.h:94
Simplex_tree()
Constructs an empty simplex tree.
Definition: Simplex_tree.h:334
void maybe_initialize_filtration()
Initializes the filtration cache if it isn't initialized yet.
Definition: Simplex_tree.h:949
boost::iterator_range< Complex_simplex_iterator > Complex_simplex_range
Range over the simplices of the simplicial complex.
Definition: Simplex_tree.h:208
Simplex_tree_boundary_simplex_iterator< Simplex_tree > Boundary_simplex_iterator
Iterator over the simplices of the boundary of a simplex.
Definition: Simplex_tree.h:194
Simplex_tree_siblings< Simplex_tree, Dictionary > Siblings
Set of nodes sharing a same parent in the simplex tree.
Definition: Simplex_tree.h:109
Simplex_handle find(const InputVertexRange &s)
Given a range of Vertex_handles, returns the Simplex_handle of the simplex in the simplicial complex ...
Definition: Simplex_tree.h:650
std::vector< Simplex_handle > Cofaces_simplex_range
Range over the cofaces of a simplex.
Definition: Simplex_tree.h:190
Boundary_opposite_vertex_simplex_range boundary_opposite_vertex_simplex_range(SimplexHandle sh)
Given a simplex, returns a range over the simplices of its boundary and their opposite vertices.
Definition: Simplex_tree.h:324
bool is_empty() const
Returns whether the complex is empty.
Definition: Simplex_tree.h:581
Simplex_handle edge_with_same_filtration(Simplex_handle sh)
Returns an edge of sh that has the same filtration value as sh if it exists, and null_simplex() other...
Definition: Simplex_tree.h:1728
Simplex_tree_complex_simplex_iterator< Simplex_tree > Complex_simplex_iterator
Iterator over the simplices of the simplicial complex.
Definition: Simplex_tree.h:206
Simplex_handle simplex(Simplex_key idx) const
Returns the simplex that has index idx in the filtration.
Definition: Simplex_tree.h:530
void initialize_filtration()
Initializes the filtration cache, i.e. sorts the simplices according to their order in the filtration...
Definition: Simplex_tree.h:925
Skeleton_simplex_range skeleton_simplex_range(int dim)
Returns a range over the simplices of the dim-skeleton of the simplicial complex.
Definition: Simplex_tree.h:255
boost::transform_iterator< return_first, Dictionary_it > Complex_vertex_iterator
Iterator over the vertices of the simplicial complex.
Definition: Simplex_tree.h:180
void insert_batch_vertices(VertexRange const &vertices, Filtration_value filt=0)
Inserts several vertices.
Definition: Simplex_tree.h:1176
std::vector< Simplex_handle > Filtration_simplex_range
Range over the simplices of the simplicial complex, ordered by the filtration.
Definition: Simplex_tree.h:218
Filtration_simplex_range::const_iterator Filtration_simplex_iterator
Iterator over the simplices of the simplicial complex, ordered by the filtration.
Definition: Simplex_tree.h:222
Extended_filtration_data extend_filtration()
Extend filtration for computing extended persistence. This function only uses the filtration values a...
Definition: Simplex_tree.h:1652
Simplex_handle minimal_simplex_with_same_filtration(Simplex_handle sh)
Returns a minimal face of sh that has the same filtration value as sh.
Definition: Simplex_tree.h:1761
Options::Vertex_handle Vertex_handle
Type for the vertex handle.
Definition: Simplex_tree.h:98
std::pair< Filtration_value, Extended_simplex_type > decode_extended_filtration(Filtration_value f, const Extended_filtration_data &efd)
Retrieve the original filtration value for a given simplex in the Simplex_tree. Since the computation...
Definition: Simplex_tree.h:1623
size_t num_vertices() const
Returns the number of vertices in the complex.
Definition: Simplex_tree.h:576
void expansion(int max_dim)
Expands the Simplex_tree containing only its one skeleton until dimension max_dim.
Definition: Simplex_tree.h:1194
Simplex_tree(const Simplex_tree &complex_source)
User-defined copy constructor reproduces the whole tree structure.
Definition: Simplex_tree.h:341
boost::iterator_range< Complex_vertex_iterator > Complex_vertex_range
Range over the vertices of the simplicial complex.
Definition: Simplex_tree.h:182
static Simplex_key null_key()
Returns a fixed number not in the interval [0, num_simplices()).
Definition: Simplex_tree.h:565
Boundary_simplex_range boundary_simplex_range(SimplexHandle sh)
Returns a range over the simplices of the boundary of a simplex.
Definition: Simplex_tree.h:307
int dimension(Simplex_handle sh)
Returns the dimension of a simplex.
Definition: Simplex_tree.h:609
bool prune_above_filtration(Filtration_value filtration)
Prune above filtration value given as parameter.
Definition: Simplex_tree.h:1464
int upper_bound_dimension() const
Returns an upper bound on the dimension of the simplicial complex.
Definition: Simplex_tree.h:620
Siblings * root()
Definition: Simplex_tree.h:902
int dimension()
Returns the dimension of the simplicial complex.
Definition: Simplex_tree.h:628
std::pair< Simplex_handle, bool > insert_simplex_and_subfaces(const InputVertexRange &Nsimplex, Filtration_value filtration=0)
Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of Vertex_handles,...
Definition: Simplex_tree.h:818
size_t num_simplices()
returns the number of simplices in the simplex_tree.
Definition: Simplex_tree.h:587
Simplex_tree_skeleton_simplex_iterator< Simplex_tree > Skeleton_simplex_iterator
Iterator over the simplices of the skeleton of the simplicial complex, for a given dimension.
Definition: Simplex_tree.h:213
Simplex_tree & operator=(Simplex_tree &&complex_source)
User-defined move assignment relocates the whole tree structure.
Definition: Simplex_tree.h:385
boost::iterator_range< Boundary_opposite_vertex_simplex_iterator > Boundary_opposite_vertex_simplex_range
Range over the simplices of the boundary of a simplex and their opposite vertices.
Definition: Simplex_tree.h:202
Simplex_tree & operator=(const Simplex_tree &complex_source)
User-defined copy assignment reproduces the whole tree structure.
Definition: Simplex_tree.h:368
void insert_graph(const OneSkeletonGraph &skel_graph)
Inserts a 1-skeleton in an empty Simplex_tree.
Definition: Simplex_tree.h:1119
boost::iterator_range< Skeleton_simplex_iterator > Skeleton_simplex_range
Range over the simplices of the skeleton of the simplicial complex, for a given dimension.
Definition: Simplex_tree.h:216
void print_hasse(std::ostream &os)
Write the hasse diagram of the simplicial complex in os.
Definition: Simplex_tree.h:1392
~Simplex_tree()
Destructor; deallocates the whole tree structure.
Definition: Simplex_tree.h:363
static Simplex_handle null_simplex()
Returns a Simplex_handle different from all Simplex_handles associated to the simplices in the simpli...
Definition: Simplex_tree.h:560
Complex_simplex_range complex_simplex_range()
Returns a range over the simplices of the simplicial complex.
Definition: Simplex_tree.h:241
Complex_vertex_range complex_vertex_range()
Returns a range over the vertices of the simplicial complex. The order is increasing according to < o...
Definition: Simplex_tree.h:230
void set_dimension(int dimension, bool exact=true)
Set a dimension for the simplicial complex.
Definition: Simplex_tree.h:912
Graph simplicial complex methods.
std::ostream & operator<<(std::ostream &os, const Permutahedral_representation< Vertex, OrderedSetPartition > &simplex)
Print a permutahedral representation to a stream.
Definition: Permutahedral_representation.h:173
This file includes common file reader for GUDHI.
bool read_simplex(std::istream &in_, std::vector< Vertex_handle > &simplex, Filtration_value &fil)
Read a face from a file.
Definition: reader_utils.h:158
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20
Node of a simplex tree with filtration value and simplex key.
Definition: Simplex_tree_node_explicit_storage.h:29
Definition: Simplex_tree.h:2012
Tag for a linear ordering of simplices.
Definition: indexing_tag.h:20
Concept describing an indexing scheme (see FilteredComplex) for applying continuous maps to a cell co...
Definition: IndexingTag.h:18
Key type used as simplex identifier.
Definition: SimplexKey.h:15
Concept of the template parameter for the class Gudhi::Simplex_tree<SimplexTreeOptions>.
Definition: SimplexTreeOptions.h:15
static const bool store_filtration
If true, each simplex has extra storage for one Filtration_value, and this value is propagated by ope...
Definition: SimplexTreeOptions.h:27
static constexpr bool contiguous_vertices
If true, the list of vertices present in the complex must always be 0, ..., num_vertices-1,...
Definition: SimplexTreeOptions.h:29
Handle type for the vertices of a cell complex.
Definition: VertexHandle.h:15