diff options
Diffstat (limited to 'src/Skeleton_blocker/include')
20 files changed, 5420 insertions, 0 deletions
diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h new file mode 100644 index 00000000..bcca851f --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker.h @@ -0,0 +1,238 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_H_ +#define SKELETON_BLOCKER_H_ + +#include <gudhi/Skeleton_blocker_complex.h> +#include <gudhi/Skeleton_blocker_geometric_complex.h> +#include <gudhi/Skeleton_blocker_simplifiable_complex.h> +#include <gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h> + +#include <gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h> +#include <gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h> + +#include <gudhi/Debug_utils.h> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** \defgroup skbl Skeleton-Blocker +@{ + +\author David Salinas + +\section skblintroduction Introduction +The Skeleton-Blocker data-structure proposes a light encoding for simplicial complexes by storing only an *implicit* representation of its +simplices +\cite socg_blockers_2011,\cite blockers2012. +Intuitively, it just stores the 1-skeleton of a simplicial complex with a graph and the set of its "missing faces" that +is very small in practice (see next section for a formal definition). +This data-structure handles all simplicial complexes operations such as + as simplex enumeration or simplex removal but operations that are particularly efficient + are operations that do not require simplex enumeration such as edge iteration, link computation or simplex contraction. + + +\section skbldefinitions Definitions + +We recall briefly classical definitions of simplicial complexes + \cite Munkres-elementsalgtop1984. +An abstract simplex is a finite non-empty set and its dimension is its number of elements minus 1. +Whenever \f$\tau \subset \sigma\f$ and \f$\tau \neq \emptyset \f$, \f$ \tau \f$ is called a face of +\f$ \sigma\f$ and \f$ \sigma\f$ is called a coface of \f$ \tau \f$ . Furthermore, +when \f$ \tau \neq \sigma\f$ we say that \f$ \tau\f$ is a proper-face of \f$ \sigma\f$. +An abstract simplicial complex is a set of simplices that contains all the faces of its simplices. +The 1-skeleton of a simplicial complex (or its graph) consists of its elements of dimension lower than 2. + + *\image html "ds_representation.png" "Skeleton-blocker representation" width=20cm + + +To encode, a simplicial complex, one can encodes all its simplices. +In case when this number gets too large, +a lighter and implicit version consists of encoding only its graph plus some elements called missing faces or blockers. +A blocker is a simplex of dimension greater than 1 +that does not belong to the complex but whose all proper faces does. + + +Remark that for a clique complex (i.e. a simplicial complex whose simplices are cliques of its graph), the set of blockers +is empty and the data-structure is then particularly sparse. +One famous example of clique-complex is the Rips complex which is intensively used +in topological data-analysis. +In practice, the set of blockers of a simplicial complex +remains also small when simplifying a Rips complex with edge contractions +but also for most of the simplicial complexes used in topological data-analysis such as Delaunay, Cech or Witness complexes. +For instance, the numbers of blockers is depicted for random 3-dimensional spheres embedded into \f$R^4\f$ +in next figure. Storing the graph and blockers of such simplicial complexes is much compact in this case than storing +their simplices. + + + *\image html "blockers_curve.png" "Number of blockers of random triangulations of 3-spheres" width=10cm + + + + +\section API + +\subsection Overview + +Two main classes of this package are Skeleton_blocker_complex and Skeleton_blocker_geometric_complex. +The first one can be used to represent an abstract simplicial complex and supports most used +operations in a simplicial complex such as : + +\li vertex/edge/simplex enumeration +\li simplifications operations such as remove star, add star (e.g. general form of collapse), +edge contractions + +The class Skeleton_blocker_geometric_complex supports the same methods as Skeleton_blocker_complex +and point access in addition. + + + +\subsection skblvisitor Visitor + +The class Skeleton_blocker_complex has a visitor that is called when usual operations such as adding an edge or remove a vertex are called. +You may want to use this visitor to compute statistics or to update another data-structure (for instance this visitor is heavily used in the \ref contr package). + + + + +\section skblexample Example + + +\subsection Iterating Iterating through vertices, edges, blockers and simplices + +Iteration through vertices, edges, simplices or blockers is straightforward with c++11 for range loops. +Note that simplex iteration with this implicit data-structure just takes +a few more time compared to iteration via an explicit representation +such as the Simplex Tree. The following example computes the Euler Characteristic +of a simplicial complex. + + \code{.cpp} + typedef Skeleton_blocker_complex<Skeleton_blocker_simple_traits> Complex; + typedef Complex::Vertex_handle Vertex_handle; + typedef Complex::Simplex Simplex; + + const int n = 15; + + // build a full complex with 10 vertices and 2^n-1 simplices + Complex complex; + for(int i=0;i<n;i++) + complex.add_vertex(); + for(int i=0;i<n;i++) + for(int j=0;j<i;j++) + complex.add_edge_without_blockers(Vertex_handle(i),Vertex_handle(j)); + + // this is just to illustrate iterators, to count number of vertices + // or edges, complex.num_vertices() and complex.num_edges() are + // more appropriated! + unsigned num_vertices = 0; + for(auto v : complex.vertex_range()){ + ++num_vertices; + } + + unsigned num_edges = 0; + for(auto e : complex.edge_range()) + ++num_edges; + + unsigned euler = 0; + unsigned num_simplices = 0; + // we use a reference to a simplex instead of a copy + // value here because a simplex is a set of integers + // and copying it cost time + for(const Simplex & s : complex.star_simplex_range()){ + ++num_simplices; + if(s.dimension()%2 == 0) + euler += 1; + else + euler -= 1; + } + std::cout << "Saw "<<num_vertices<<" vertices, "<<num_edges<<" edges and "<<num_simplices<<" simplices"<<std::endl; + std::cout << "The Euler Characteristic is "<<euler<<std::endl; + \endcode + + +\verbatim +./SkeletonBlockerIteration +Saw 15 vertices, 105 edges and 32767 simplices +The Euler Characteristic is 1 + 0.537302s wall, 0.530000s user + 0.000000s system = 0.530000s CPU (98.6%) +\endverbatim + + +\subsection s Constructing a skeleton-blockers from a list of maximal faces or from a list of faces + + \code{.cpp} + std::vector<Simplex> simplices; + + //add 4 triangles of a tetrahedron 0123 + simplices.push_back(Simplex(Vertex_handle(0),Vertex_handle(1),Vertex_handle(2))); + simplices.push_back(Simplex(Vertex_handle(1),Vertex_handle(2),Vertex_handle(3))); + simplices.push_back(Simplex(Vertex_handle(3),Vertex_handle(0),Vertex_handle(2))); + simplices.push_back(Simplex(Vertex_handle(3),Vertex_handle(0),Vertex_handle(1))); + + Complex complex; + //get complex from top faces + make_complex_from_top_faces(complex,simplices.begin(),simplices.end()); + + std::cout << "Simplices:"<<std::endl; + for(const Simplex & s : complex.star_simplex_range()) + std::cout << s << " "; + std::cout << std::endl; + + //One blocker as simplex 0123 is not in the complex but all its proper faces are. + std::cout << "Blockers: "<<complex.blockers_to_string()<<std::endl; + + //now build a complex from its full list of simplices + simplices.clear(); + simplices.push_back(Simplex(Vertex_handle(0))); + simplices.push_back(Simplex(Vertex_handle(1))); + simplices.push_back(Simplex(Vertex_handle(2))); + simplices.push_back(Simplex(Vertex_handle(0),Vertex_handle(1))); + simplices.push_back(Simplex(Vertex_handle(1),Vertex_handle(2))); + simplices.push_back(Simplex(Vertex_handle(2),Vertex_handle(0))); + complex = Complex(simplices.begin(),simplices.end()); + + std::cout << "Simplices:"<<std::endl; + for(const Simplex & s : complex.star_simplex_range()) + std::cout << s << " "; + std::cout << std::endl; + + //One blocker as simplex 012 is not in the complex but all its proper faces are. + std::cout << "Blockers: "<<complex.blockers_to_string()<<std::endl; + \endcode +\verbatim +./SkeletonBlockerFromSimplices +Simplices: +{0} {0,1} {0,2} {0,3} {0,1,2} {0,1,3} {0,2,3} {1} {1,2} {1,3} {1,2,3} {2} {2,3} {3} +Blockers: {0,1,2,3} + +Simplices: +{0} {0,1} {0,2} {1} {1,2} {2} +Blockers: {0,1,2} +\endverbatim + + +\section Acknowledgements +The author wishes to thank Dominique Attali and André Lieutier for +their collaboration to write the two initial papers +\cite socg_blockers_2011,\cite blockers2012 + about this data-structure + and also Dominique for leaving him use a prototype. + +@} */ + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h new file mode 100644 index 00000000..9f145013 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h @@ -0,0 +1,132 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_SKELETON_BLOCKER_COMPLEX_VISITOR_H_ +#define SKELETON_BLOCKER_SKELETON_BLOCKER_COMPLEX_VISITOR_H_ + +#include <gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h> + +namespace Gudhi { + +namespace skeleton_blocker { +// TODO(DS): to be constified + +/** + *@class Skeleton_blocker_complex_visitor + *@brief Interface for a visitor of a simplicial complex. + */ +template<typename Vertex_handle> +class Skeleton_blocker_complex_visitor { + public: + virtual ~Skeleton_blocker_complex_visitor() { } + + virtual void on_add_vertex(Vertex_handle) = 0; + virtual void on_remove_vertex(Vertex_handle) = 0; + + virtual void on_add_edge_without_blockers(Vertex_handle a, Vertex_handle b) = 0; + virtual void on_remove_edge(Vertex_handle a, Vertex_handle b) = 0; + + /** + * @brief Called when an edge changes, during the contraction of + * an edge + */ + virtual void on_changed_edge(Vertex_handle a, Vertex_handle b) = 0; + + /** + * @brief Called when performing an edge contraction when + * an edge bx is replaced by an edge ax (not already present). + * Precisely, this methods is called this way in contract_edge : + * add_edge_without_blockers(a,x) + * on_swaped_edge(a,b,x) + * remove_edge(b,x) + */ + virtual void on_swaped_edge(Vertex_handle a, Vertex_handle b, + Vertex_handle x) = 0; + virtual void on_add_blocker( + const Skeleton_blocker_simplex<Vertex_handle>&) = 0; + virtual void on_delete_blocker( + const Skeleton_blocker_simplex<Vertex_handle>*) = 0; +}; + +/** + *@class Dummy_complex_visitor + *@brief A dummy visitor of a simplicial complex that does nothing + * + */ +template<typename Vertex_handle> +class Dummy_complex_visitor : public Skeleton_blocker_complex_visitor< +Vertex_handle> { + public: + void on_add_vertex(Vertex_handle) { } + + void on_remove_vertex(Vertex_handle) { } + + void on_add_edge_without_blockers(Vertex_handle a, Vertex_handle b) { } + + void on_remove_edge(Vertex_handle a, Vertex_handle b) { } + + void on_changed_edge(Vertex_handle a, Vertex_handle b) { } + + void on_swaped_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x) { } + + void on_add_blocker(const Skeleton_blocker_simplex<Vertex_handle>&) { } + + void on_delete_blocker(const Skeleton_blocker_simplex<Vertex_handle>*) { } +}; + +/** + *@class Print_complex_visitor + *@brief A visitor that prints the visit information. + * + */ +template<typename Vertex_handle> +class Print_complex_visitor : public Skeleton_blocker_complex_visitor< +Vertex_handle> { + public: + void on_add_vertex(Vertex_handle v) { + std::cerr << "on_add_vertex:" << v << std::endl; + } + + void on_remove_vertex(Vertex_handle v) { + std::cerr << "on_remove_vertex:" << v << std::endl; + } + + void on_add_edge_without_blockers(Vertex_handle a, Vertex_handle b) { + std::cerr << "on_add_edge_without_blockers:" << a << "," << b << std::endl; + } + + void on_remove_edge(Vertex_handle a, Vertex_handle b) { + std::cerr << "on_remove_edge:" << a << "," << b << std::endl; + } + + void on_changed_edge(Vertex_handle a, Vertex_handle b) { + std::cerr << "on_changed_edge:" << a << "," << b << std::endl; + } + + void on_swaped_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x) { + std::cerr << "on_swaped_edge:" << a << "," << b << "," << x << std::endl; + } + + void on_add_blocker(const Skeleton_blocker_simplex<Vertex_handle>& b) { + std::cerr << "on_add_blocker:" << b << std::endl; + } + + void on_delete_blocker(const Skeleton_blocker_simplex<Vertex_handle>* b) { + std::cerr << "on_delete_blocker:" << b << std::endl; + } +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_SKELETON_BLOCKER_COMPLEX_VISITOR_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h new file mode 100644 index 00000000..d348b696 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h @@ -0,0 +1,65 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_SKELETON_BLOCKER_LINK_SUPERIOR_H_ +#define SKELETON_BLOCKER_SKELETON_BLOCKER_LINK_SUPERIOR_H_ + +#include <gudhi/Skeleton_blocker_link_complex.h> + +namespace Gudhi { + +namespace skeleton_blocker { + +template<class ComplexType> class Skeleton_blocker_sub_complex; + +/** + * \brief Class representing the link of a simplicial complex encoded by a skeleton/blockers pair. + * It computes only vertices greater than the simplex used to build the link. + */ +template<typename ComplexType> +class Skeleton_blocker_link_superior : public Skeleton_blocker_link_complex< +ComplexType> { + typedef typename ComplexType::Edge_handle Edge_handle; + + typedef typename ComplexType::boost_vertex_handle boost_vertex_handle; + + public: + typedef typename ComplexType::Vertex_handle Vertex_handle; + typedef typename ComplexType::Root_vertex_handle Root_vertex_handle; + typedef typename ComplexType::Simplex Simplex; + typedef typename ComplexType::Root_simplex_handle Root_simplex_handle; + typedef typename ComplexType::BlockerMap BlockerMap; + typedef typename ComplexType::BlockerPair BlockerPair; + typedef typename ComplexType::BlockerMapIterator BlockerMapIterator; + typedef typename ComplexType::BlockerMapConstIterator BlockerMapConstIterator; + typedef typename ComplexType::Simplex::Simplex_vertex_const_iterator AddressSimplexConstIterator; + typedef typename ComplexType::Root_simplex_handle::Simplex_vertex_const_iterator IdSimplexConstIterator; + + Skeleton_blocker_link_superior() + : Skeleton_blocker_link_complex<ComplexType>(true) { } + + Skeleton_blocker_link_superior(const ComplexType & parent_complex, + Simplex& alpha_parent_adress) + : Skeleton_blocker_link_complex<ComplexType>(parent_complex, + alpha_parent_adress, true) { } + + Skeleton_blocker_link_superior(const ComplexType & parent_complex, + Vertex_handle a_parent_adress) + : Skeleton_blocker_link_complex<ComplexType>(parent_complex, + a_parent_adress, true) { } +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_SKELETON_BLOCKER_LINK_SUPERIOR_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h new file mode 100644 index 00000000..52300493 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_off_io.h @@ -0,0 +1,191 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_SKELETON_BLOCKER_OFF_IO_H_ +#define SKELETON_BLOCKER_SKELETON_BLOCKER_OFF_IO_H_ + +#include <gudhi/Off_reader.h> + +#include <string> +#include <vector> +#include <map> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + *@brief Off reader visitor that can be passed to Off_reader to read a Skeleton_blocker_complex. + */ +template<typename Complex> +class Skeleton_blocker_off_flag_visitor_reader { + Complex& complex_; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Point Point; + + const bool load_only_points_; + + public: + explicit Skeleton_blocker_off_flag_visitor_reader(Complex& complex, bool load_only_points = false) : + complex_(complex), + load_only_points_(load_only_points) { } + + void init(int dim, int num_vertices, int num_faces, int num_edges) { + // TODO(DS): do an assert to check that this number are correctly read + // TODO(DS): reserve size for vector points + } + + void point(const std::vector<double>& point) { + complex_.add_vertex(Point(point.begin(), point.end())); + } + + void maximal_face(const std::vector<int>& face) { + if (!load_only_points_) { + for (size_t i = 0; i < face.size(); ++i) + for (size_t j = i + 1; j < face.size(); ++j) { + complex_.add_edge_without_blockers(Vertex_handle(face[i]), Vertex_handle(face[j])); + } + } + } + + void done() { } +}; + +/** + *@brief Off reader visitor that can be passed to Off_reader to read a Skeleton_blocker_complex. + */ +template<typename Complex> +class Skeleton_blocker_off_visitor_reader { + Complex& complex_; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Simplex Simplex; + typedef typename Complex::Point Point; + + const bool load_only_points_; + std::vector<Point> points_; + std::vector<Simplex> maximal_faces_; + + public: + explicit Skeleton_blocker_off_visitor_reader(Complex& complex, bool load_only_points = false) : + complex_(complex), + load_only_points_(load_only_points) { } + + void init(int dim, int num_vertices, int num_faces, int num_edges) { + maximal_faces_.reserve(num_faces); + points_.reserve(num_vertices); + } + + void point(const std::vector<double>& point) { + points_.emplace_back(point.begin(), point.end()); + } + + void maximal_face(const std::vector<int>& face) { + if (!load_only_points_) { + Simplex s; + for (auto x : face) + s.add_vertex(Vertex_handle(x)); + maximal_faces_.emplace_back(s); + } + } + + void done() { + complex_ = make_complex_from_top_faces<Complex>(maximal_faces_.begin(), maximal_faces_.end(), + points_.begin(), points_.end()); + } +}; + +/** + *@brief Class that allows to load a Skeleton_blocker_complex from an off file. + */ +template<typename Complex> +class Skeleton_blocker_off_reader { + public: + /** + * name_file : file to read + * read_complex : complex that will receive the file content + * read_only_points : specify true if only the points must be read + */ + Skeleton_blocker_off_reader(const std::string & name_file, Complex& read_complex, + bool read_only_points = false, bool is_flag = false) : valid_(false) { + std::ifstream stream(name_file); + if (stream.is_open()) { + if (is_flag || read_only_points) { + Skeleton_blocker_off_flag_visitor_reader<Complex> off_visitor(read_complex, read_only_points); + Off_reader off_reader(stream); + valid_ = off_reader.read(off_visitor); + } else { + Skeleton_blocker_off_visitor_reader<Complex> off_visitor(read_complex, read_only_points); + Off_reader off_reader(stream); + valid_ = off_reader.read(off_visitor); + } + } + } + + /** + * return true if reading did not meet problems. + */ + bool is_valid() const { + return valid_; + } + + private: + bool valid_; +}; + +template<typename Complex> +class Skeleton_blocker_off_writer { + public: + /** + * name_file : file where the off will be written + * save_complex : complex that be outputted in the file + * for now only save triangles. + */ + Skeleton_blocker_off_writer(const std::string & name_file, const Complex& save_complex) { + typedef typename Complex::Vertex_handle Vertex_handle; + + std::ofstream stream(name_file); + if (stream.is_open()) { + stream << "OFF\n"; + size_t num_triangles = std::distance(save_complex.triangle_range().begin(), save_complex.triangle_range().end()); + stream << save_complex.num_vertices() << " " << num_triangles << " 0 \n"; + + // in case the complex has deactivated some vertices, eg only has vertices 0 2 5 7 for instance + // we compute a map from 0 2 5 7 to 0 1 2 3 + std::map<Vertex_handle, size_t> vertex_num; + size_t current_vertex = 0; + + for (auto v : save_complex.vertex_range()) { + vertex_num[v] = current_vertex++; + const auto& pt(save_complex.point(v)); + for (auto x : pt) + stream << x << " "; + stream << std::endl; + } + + for (const auto & t : save_complex.triangle_range()) { + stream << "3 "; + for (auto x : t) + stream << vertex_num[x] << " "; + stream << std::endl; + } + stream.close(); + } else { + std::cerr << "could not open file " << name_file << std::endl; + } + } +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_SKELETON_BLOCKER_OFF_IO_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h new file mode 100644 index 00000000..772e33aa --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h @@ -0,0 +1,87 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_SKELETON_BLOCKER_SIMPLE_GEOMETRIC_TRAITS_H_ +#define SKELETON_BLOCKER_SKELETON_BLOCKER_SIMPLE_GEOMETRIC_TRAITS_H_ + +#include <gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h> + +#include <string> +#include <sstream> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + * @extends SkeletonBlockerGeometricDS + * @ingroup skbl + * @brief Simple traits that is a model of SkeletonBlockerGeometricDS and + * can be passed as a template argument to Skeleton_blocker_geometric_complex + */ +template<typename GeometryTrait> +struct Skeleton_blocker_simple_geometric_traits : +public Skeleton_blocker_simple_traits { + public: + typedef GeometryTrait GT; + typedef typename GT::Point Point; + typedef typename Skeleton_blocker_simple_traits::Root_vertex_handle Root_vertex_handle; + typedef typename Skeleton_blocker_simple_traits::Graph_vertex Simple_vertex; + + /** + * @brief Vertex with a point attached. + */ + class Simple_geometric_vertex : public Simple_vertex { + template<class ComplexGeometricTraits> friend class Skeleton_blocker_geometric_complex; + private: + Point point_; + + Point& point() { + return point_; + } + + const Point& point() const { + return point_; + } + }; + + class Simple_geometric_edge : + public Skeleton_blocker_simple_traits::Graph_edge { + int index_; + public: + Simple_geometric_edge() + : Skeleton_blocker_simple_traits::Graph_edge(), + index_(-1) { } + + /** + * @brief Allows to modify the index of the edge. + * The indices of the edge are used to store heap information + * in the edge contraction algorithm. + */ + int& index() { + return index_; + } + + int index() const { + return index_; + } + }; + + typedef Simple_geometric_vertex Graph_vertex; + typedef Skeleton_blocker_simple_traits::Graph_edge Graph_edge; +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_SKELETON_BLOCKER_SIMPLE_GEOMETRIC_TRAITS_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h new file mode 100644 index 00000000..0c0cc624 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simple_traits.h @@ -0,0 +1,178 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_SKELETON_BLOCKER_SIMPLE_TRAITS_H_ +#define SKELETON_BLOCKER_SKELETON_BLOCKER_SIMPLE_TRAITS_H_ + +#include <gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h> + +#include <string> +#include <sstream> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + * @extends SkeletonBlockerDS + * @ingroup skbl + * @brief Simple traits that is a model of SkeletonBlockerDS and + * can be passed as a template argument to Skeleton_blocker_complex + */ +struct Skeleton_blocker_simple_traits { + /** + * @brief Global and local handle similar to <a href="http://www.boost.org/doc/libs/1_38_0/libs/graph/doc/subgraph.html">boost subgraphs</a>. + * Vertices are stored in a vector. + * For the root simplicial complex, the local and global descriptors are the same. + * For a subcomplex L and one of its vertices 'v', the local descriptor of 'v' is its position in + * the vertex vector of the subcomplex L whereas its global descriptor is the position of 'v' + * in the vertex vector of the root simplicial complex. + */ + struct Root_vertex_handle { + typedef int boost_vertex_handle; + + explicit Root_vertex_handle(boost_vertex_handle val = -1) + : vertex(val) { } + boost_vertex_handle vertex; + + bool operator!=(const Root_vertex_handle& other) const { + return !(this->vertex == other.vertex); + } + + bool operator==(const Root_vertex_handle& other) const { + return this->vertex == other.vertex; + } + + bool operator<(const Root_vertex_handle& other) const { + return this->vertex < other.vertex; + } + + friend std::ostream& operator<<(std::ostream& o, + const Root_vertex_handle & v) { + o << v.vertex; + return o; + } + }; + + struct Vertex_handle { + typedef int boost_vertex_handle; + + explicit Vertex_handle(boost_vertex_handle val = -1) + : vertex(val) { } + + operator int() const { + return static_cast<int> (vertex); + } + + boost_vertex_handle vertex; + + bool operator==(const Vertex_handle& other) const { + return this->vertex == other.vertex; + } + + bool operator!=(const Vertex_handle& other) const { + return this->vertex != other.vertex; + } + + bool operator<(const Vertex_handle& other) const { + return this->vertex < other.vertex; + } + + friend std::ostream& operator<<(std::ostream& o, const Vertex_handle & v) { + o << v.vertex; + return o; + } + }; + + class Graph_vertex { + bool is_active_; + Root_vertex_handle id_; + + public: + virtual ~Graph_vertex() { } + + void activate() { + is_active_ = true; + } + + void deactivate() { + is_active_ = false; + } + + bool is_active() const { + return is_active_; + } + + void set_id(Root_vertex_handle i) { + id_ = i; + } + + Root_vertex_handle get_id() const { + return id_; + } + + virtual std::string to_string() const { + std::ostringstream res; + res << id_; + return res.str(); + } + + friend std::ostream& operator<<(std::ostream& o, const Graph_vertex & v) { + o << v.to_string(); + return o; + } + }; + + class Graph_edge { + Root_vertex_handle a_; + Root_vertex_handle b_; + int index_; + + public: + Graph_edge() + : a_(-1), + b_(-1), + index_(-1) { } + + int& index() { + return index_; + } + + int index() const { + return index_; + } + + void setId(Root_vertex_handle a, Root_vertex_handle b) { + a_ = a; + b_ = b; + } + + Root_vertex_handle first() const { + return a_; + } + + Root_vertex_handle second() const { + return b_; + } + + friend std::ostream& operator<<(std::ostream& o, const Graph_edge & v) { + o << "(" << v.a_ << "," << v.b_ << " - id = " << v.index(); + return o; + } + }; +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_SKELETON_BLOCKER_SIMPLE_TRAITS_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h new file mode 100644 index 00000000..12fe6469 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h @@ -0,0 +1,362 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_SKELETON_BLOCKER_SIMPLEX_H_ +#define SKELETON_BLOCKER_SKELETON_BLOCKER_SIMPLEX_H_ + +#include <cassert> +#include <iostream> +#include <set> +#include <vector> +#include <initializer_list> +#include <string> +#include <algorithm> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + *@brief Abstract simplex used in Skeleton blockers data-structure. + * + * An abstract simplex is represented as an ordered set of T elements, + * each element representing a vertex. + * + * The element representing a vertex can be SkeletonBlockerDS::Vertex_handle or SkeletonBlockerDS::Root_vertex_handle. + * + * + */ +template<typename T> + +class Skeleton_blocker_simplex { + private: + std::set<T> simplex_set; + + public: + typedef typename T::boost_vertex_handle boost_vertex_handle; + + typedef T Vertex_handle; + + typedef typename std::set<T>::const_iterator Simplex_vertex_const_iterator; + typedef typename std::set<T>::iterator Simplex_vertex_iterator; + + /** @name Constructors / Destructors / Initialization + */ + //@{ + + void clear() { + simplex_set.clear(); + } + + Skeleton_blocker_simplex(std::initializer_list<T>& list) { + std::for_each(list.begin(), list.end(), [&] (T const& elt) { + add_vertex(elt); + }); + } + + template<typename ... Args> + explicit Skeleton_blocker_simplex(Args ... args) { + add_vertices(args...); + } + + template<typename ... Args> + void add_vertices(T v, Args ... args) { + add_vertex(v); + add_vertices(args...); + } + + void add_vertices(T v) { + add_vertex(v); + } + + void add_vertices() { } + + /** + * Initialize a simplex with a string such as {0,1,2} + */ + explicit Skeleton_blocker_simplex(std::string token) { + clear(); + if ((token[0] == '{') && (token[token.size() - 1] == '}')) { + token.erase(0, 1); + token.erase(token.size() - 1, 1); + while (token.size() != 0) { + int coma_position = token.find_first_of(','); + // cout << "coma_position:"<<coma_position<<endl; + std::string n = token.substr(0, coma_position); + // cout << "token:"<<token<<endl; + token.erase(0, n.size() + 1); + add_vertex((T) (atoi(n.c_str()))); + } + } + } + + //@} + + /** @name Simplex manipulation + */ + //@{ + /** + * Add the vertex v to the simplex: + * Note that adding two times the same vertex is + * the same that adding it once. + */ + void add_vertex(T v) { + simplex_set.insert(v); + } + + /** + * Remove the vertex v from the simplex: + */ + void remove_vertex(T v) { + simplex_set.erase(v); + } + + /** + * Intersects the simplex with the simplex a. + */ + void intersection(const Skeleton_blocker_simplex & a) { + std::vector<T> v; + v.reserve((std::min)(simplex_set.size(), a.simplex_set.size())); + + set_intersection(simplex_set.begin(), simplex_set.end(), + a.simplex_set.begin(), a.simplex_set.end(), + std::back_inserter(v)); + clear(); + for (auto i : v) + simplex_set.insert(i); + } + + /** + * Substracts a from the simplex. + */ + void difference(const Skeleton_blocker_simplex & a) { + std::vector<T> v; + v.reserve(simplex_set.size()); + + set_difference(simplex_set.begin(), simplex_set.end(), + a.simplex_set.begin(), a.simplex_set.end(), + std::back_inserter(v)); + + clear(); + for (auto i : v) + simplex_set.insert(i); + } + + /** + * Add vertices of a to the simplex. + */ + void union_vertices(const Skeleton_blocker_simplex & a) { + std::vector<T> v; + v.reserve(simplex_set.size() + a.simplex_set.size()); + + set_union(simplex_set.begin(), simplex_set.end(), a.simplex_set.begin(), + a.simplex_set.end(), std::back_inserter(v)); + clear(); + simplex_set.insert(v.begin(), v.end()); + } + + typename std::set<T>::const_iterator begin() const { + return simplex_set.cbegin(); + } + + typename std::set<T>::const_iterator end() const { + return simplex_set.cend(); + } + + typename std::set<T>::const_reverse_iterator rbegin() const { + return simplex_set.crbegin(); + } + + typename std::set<T>::const_reverse_iterator rend() const { + return simplex_set.crend(); + } + + typename std::set<T>::iterator begin() { + return simplex_set.begin(); + } + + typename std::set<T>::iterator end() { + return simplex_set.end(); + } + + //@} + + /** @name Queries + */ + //@{ + /** + * Returns the dimension of the simplex. + */ + int dimension() const { + return (simplex_set.size() - 1); + } + + bool empty() const { + return simplex_set.empty(); + } + + /** + * Returns the first and smallest vertex of the simplex. + * + * Be careful : assumes the simplex is non-empty. + */ + T first_vertex() const { + assert(!empty()); + return *(begin()); + } + + /** + * Returns the last and greatest vertex of the simplex. + * + * Be careful : assumes the simplex is non-empty. + */ + T last_vertex() const { + assert(!empty()); + return *(simplex_set.rbegin()); + } + + /** + * @return true iff the simplex contains the simplex a. + */ + bool contains(const Skeleton_blocker_simplex & a) const { + return includes(simplex_set.cbegin(), simplex_set.cend(), + a.simplex_set.cbegin(), a.simplex_set.cend()); + } + + /** + * @return true iff the simplex contains the difference \f$ a \setminus b \f$. + */ + bool contains_difference(const Skeleton_blocker_simplex& a, + const Skeleton_blocker_simplex& b) const { + auto first1 = begin(); + auto last1 = end(); + + auto first2 = a.begin(); + auto last2 = a.end(); + + while (first2 != last2) { + // we ignore vertices of b + if (b.contains(*first2)) { + ++first2; + } else { + if ((first1 == last1) || (*first2 < *first1)) + return false; + if (!(*first1 < *first2)) + ++first2; + ++first1; + } + } + return true; + } + + /** + * @return true iff the simplex contains the difference \f$ a \setminus \{ x \} \f$. + */ + bool contains_difference(const Skeleton_blocker_simplex& a, T x) const { + auto first1 = begin(); + auto last1 = end(); + + auto first2 = a.begin(); + auto last2 = a.end(); + + while (first2 != last2) { + // we ignore vertices x + if (x == *first2) { + ++first2; + } else { + if ((first1 == last1) || (*first2 < *first1)) + return false; + if (!(*first1 < *first2)) + ++first2; + ++first1; + } + } + return true; + } + + /** + * @return true iff the simplex contains the difference \f$ a \setminus \{ x,y \} \f$. + */ + bool contains_difference(const Skeleton_blocker_simplex& a, T x, T y) const { + auto first1 = begin(); + auto last1 = end(); + + auto first2 = a.begin(); + auto last2 = a.end(); + + while (first2 != last2) { + // we ignore vertices of x,y + if (x == *first2 || y == *first2) { + ++first2; + } else { + if ((first1 == last1) || (*first2 < *first1)) + return false; + if (!(*first1 < *first2)) + ++first2; + ++first1; + } + } + return true; + } + + bool contains(T v) const { + return (simplex_set.find(v) != simplex_set.end()); + } + + bool disjoint(const Skeleton_blocker_simplex& a) const { + std::vector<T> v; + v.reserve(std::min(simplex_set.size(), a.simplex_set.size())); + + set_intersection(simplex_set.cbegin(), simplex_set.cend(), + a.simplex_set.cbegin(), a.simplex_set.cend(), + std::back_inserter(v)); + + return (v.size() == 0); + } + + bool operator==(const Skeleton_blocker_simplex& other) const { + return (this->simplex_set == other.simplex_set); + } + + bool operator!=(const Skeleton_blocker_simplex& other) const { + return (this->simplex_set != other.simplex_set); + } + + bool operator<(const Skeleton_blocker_simplex& other) const { + return (std::lexicographical_compare(this->simplex_set.begin(), + this->simplex_set.end(), other.begin(), + other.end())); + } + + //@} + + friend std::ostream& operator<<(std::ostream& o, + const Skeleton_blocker_simplex & sigma) { + bool first = true; + o << "{"; + for (auto i : sigma) { + if (first) + first = false; + else + o << ","; + o << i; + } + o << "}"; + return o; + } +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_SKELETON_BLOCKER_SIMPLEX_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h new file mode 100644 index 00000000..4c48ff31 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h @@ -0,0 +1,261 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_SKELETON_BLOCKER_SUB_COMPLEX_H_ +#define SKELETON_BLOCKER_SKELETON_BLOCKER_SUB_COMPLEX_H_ + +#include <gudhi/Skeleton_blocker_complex.h> +#include <gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h> +#include <gudhi/Debug_utils.h> + +#include <map> +#include <vector> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + * @brief Simplicial subcomplex of a complex represented by a skeleton/blockers pair. + * @extends Skeleton_blocker_complex + * @details Stores a subcomplex of a simplicial complex. + * To simplify explanations below, we will suppose that : + * - K is the root simplicial complex + * - L is a subcomplex of K. + * + * One vertex of K may exists in L but with a different address. + * To be able to locate the vertices in K from vertices of L, the class + * stores a map 'adresses' between vertices of K and vertices of L. + * + * Note that the type for handle of vertices of L is 'Vertex_handle' and + * the type for handle of vertices of K is 'Root_vertex_handle'. + * + * The template ComplexType is type of the root complex. It allows to know + * if the subcomplex is geometric or not. + * It has to be either 'Skeleton_blockers_complex' or 'Skeleton_blockers_geometric_complex'. + * + */ +template<typename ComplexType> +class Skeleton_blocker_sub_complex : public ComplexType { + protected: + template<class T> friend class Skeleton_blocker_link_complex; + + typedef typename ComplexType::Graph Graph; + typedef typename ComplexType::Edge_handle Edge_handle; + + typedef typename ComplexType::boost_vertex_handle boost_vertex_handle; + + public: + using ComplexType::add_vertex; + using ComplexType::add_edge_without_blockers; + using ComplexType::add_blocker; + + typedef typename ComplexType::Vertex_handle Vertex_handle; + typedef typename ComplexType::Root_vertex_handle Root_vertex_handle; + typedef typename ComplexType::Simplex Simplex; + typedef typename ComplexType::Root_simplex_handle Root_simplex_handle; + + protected: + /** + * @brief Determines whether all proper faces of simplex 'sigma' belong to 'link1' \cup 'link2' + * where 'link1' and 'link2' are subcomplexes of the same complex of type ComplexType + */ + typedef std::map<Root_vertex_handle, Vertex_handle> IdAddressMap; + typedef typename IdAddressMap::value_type AddressPair; + typedef typename IdAddressMap::iterator IdAddressMapIterator; + typedef typename IdAddressMap::const_iterator IdAddressMapConstIterator; + std::map<Root_vertex_handle, Vertex_handle> adresses; + + public: + /** + * Add a vertex 'global' of K to L. When added to L, this vertex will receive + * another number, addresses(global), its local adress. + * return the adress where the vertex lay on L. + * The vertex corresponding to 'global' must not be already present + * in the complex. + */ + Vertex_handle add_vertex(Root_vertex_handle global) { + assert(!this->contains_vertex(global)); + Vertex_handle address(boost::add_vertex(this->skeleton)); + this->num_vertices_++; + (*this)[address].activate(); + (*this)[address].set_id(global); + adresses.insert(AddressPair(global, address)); + this->degree_.push_back(0); + return address; + } + + /** + * Add an edge (v1_root,v2_root) to the sub-complex. + * It assumes that both vertices corresponding to v1_root and v2_root are present + * in the sub-complex. + */ + void add_edge_without_blockers(Root_vertex_handle v1_root, Root_vertex_handle v2_root) { + auto v1_sub(this->get_address(v1_root)); + auto v2_sub(this->get_address(v2_root)); + assert(v1_sub && v2_sub); + this->ComplexType::add_edge_without_blockers(*v1_sub, *v2_sub); + } + + /** + * Add a blocker to the sub-complex. + * It assumes that all vertices of blocker_root are present + * in the sub-complex. + */ + void add_blocker(const Root_simplex_handle& blocker_root) { + auto blocker_sub = this->get_address(blocker_root); + assert(blocker_sub); + this->add_blocker(new Simplex(*blocker_sub)); + } + + public: + /** + * Constructs the restricted complex of 'parent_complex' to + * vertices of 'simplex'. + */ + void make_restricted_complex(const ComplexType & parent_complex, + const Simplex& simplex) { + this->clear(); + // add vertices to the sub complex + for (auto x : simplex) { + assert(parent_complex.contains_vertex(x)); + auto x_local = this->add_vertex(parent_complex[x].get_id()); + (*this)[x_local] = parent_complex[x]; + } + + // add edges to the sub complex + for (auto x : simplex) { + // x_neigh is the neighbor of x intersected with vertices_simplex + Simplex x_neigh; + parent_complex.add_neighbours(x, x_neigh, true); + x_neigh.intersection(simplex); + for (auto y : x_neigh) { + this->add_edge_without_blockers(parent_complex[x].get_id(), parent_complex[y].get_id()); + } + } + + // add blockers to the sub complex + for (auto blocker : parent_complex.const_blocker_range()) { + // check if it is the first time we encounter the blocker + if (simplex.contains(*blocker)) { + Root_simplex_handle blocker_root(parent_complex.get_id(*(blocker))); + Simplex blocker_restr( + *(this->get_simplex_address(blocker_root))); + this->add_blocker(new Simplex(blocker_restr)); + } + } + } + + void clear() { + adresses.clear(); + ComplexType::clear(); + } + + /** + * Compute the local vertex in L corresponding to the vertex global in K. + * runs in O(log n) if n = num_vertices() + */ + boost::optional<Vertex_handle> get_address(Root_vertex_handle global) const { + boost::optional < Vertex_handle > res; + IdAddressMapConstIterator it = adresses.find(global); + if (it == adresses.end()) + res.reset(); + else + res = (*it).second; + return res; + } + + // /** + // * Allocates a simplex in L corresponding to the simplex s in K + // * with its local adresses and returns an AddressSimplex. + // */ + // boost::optional<Simplex> get_address(const Root_simplex_handle & s) const; + + // private: + /** + * same as get_address except that it will return a simplex in any case. + * The vertices that were not found are not added. + */ + // @remark should be private but problem with VS + + std::vector<boost::optional<Vertex_handle> > get_addresses( + const Root_simplex_handle & s) const { + std::vector < boost::optional<Vertex_handle> > res; + for (auto i : s) { + res.push_back(get_address(i)); + } + return res; + } +}; + +/** + * @remark remarque perte de temps a creer un nouveau simplexe a chaque fois + * alors qu'on pourrait utiliser a la place de 'addresses_sigma_in_link' + * un simplex avec des valeurs sp�ciales ComplexDS::null_vertex par exemple + * pour indiquer qu'un vertex n'appartient pas au complex + */ +template<typename ComplexType> +bool proper_face_in_union( + Skeleton_blocker_sub_complex<ComplexType> & link, + std::vector<boost::optional<typename ComplexType::Vertex_handle> > & addresses_sigma_in_link, + std::size_t vertex_to_be_ignored) { + // we test that all vertices of 'addresses_sigma_in_link' but 'vertex_to_be_ignored' + // are in link1 if it is the case we construct the corresponding simplex + bool vertices_sigma_are_in_link = true; + typename ComplexType::Simplex sigma_in_link; + for (std::size_t i = 0; i < addresses_sigma_in_link.size(); ++i) { + if (i != vertex_to_be_ignored) { + if (!addresses_sigma_in_link[i]) { + vertices_sigma_are_in_link = false; + break; + } else { + sigma_in_link.add_vertex(*addresses_sigma_in_link[i]); + } + } + } + // If one of vertices of the simplex is not in the complex then it returns false + // Otherwise, it tests if the simplex is in the complex + return vertices_sigma_are_in_link && link.contains(sigma_in_link); +} + +// Remark: this function should be friend in order to leave get_adresses private +// however doing so seemes currently not possible due to a visual studio bug c2668 +// "the compiler does not support partial ordering of template functions as specified in the C++ Standard" +// http://www.serkey.com/error-c2668-ambiguous-call-to-overloaded-function-bb45ft.html + +template<typename ComplexType> +bool proper_faces_in_union( + Skeleton_blocker_simplex<typename ComplexType::Root_vertex_handle> & sigma, + Skeleton_blocker_sub_complex<ComplexType> & link1, + Skeleton_blocker_sub_complex<ComplexType> & link2) { + typedef typename ComplexType::Vertex_handle Vertex_handle; + std::vector < boost::optional<Vertex_handle> > addresses_sigma_in_link1 = + link1.get_addresses(sigma); + std::vector < boost::optional<Vertex_handle> > addresses_sigma_in_link2 = + link2.get_addresses(sigma); + + for (std::size_t current_index = 0; current_index < addresses_sigma_in_link1.size(); + ++current_index) { + if (!proper_face_in_union(link1, addresses_sigma_in_link1, current_index) + && !proper_face_in_union(link2, addresses_sigma_in_link2, + current_index)) { + return false; + } + } + return true; +} + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_SKELETON_BLOCKER_SUB_COMPLEX_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Top_faces.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Top_faces.h new file mode 100644 index 00000000..91e79b42 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Top_faces.h @@ -0,0 +1,61 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_INTERNAL_TOP_FACES_H_ +#define SKELETON_BLOCKER_INTERNAL_TOP_FACES_H_ + +#include <list> +#include <vector> +#include <set> + +namespace Gudhi { + +namespace skeleton_blocker { + +template<typename SimplexHandle> +std::list<SimplexHandle> subfaces(SimplexHandle top_face) { + std::list<SimplexHandle> res; + if (top_face.dimension() == -1) return res; + if (top_face.dimension() == 0) { + res.push_back(top_face); + return res; + } else { + auto first_vertex = top_face.first_vertex(); + top_face.remove_vertex(first_vertex); + res = subfaces(top_face); + std::list<SimplexHandle> copy = res; + for (auto& simplex : copy) { + simplex.add_vertex(first_vertex); + } + res.push_back(SimplexHandle(first_vertex)); + res.splice(res.end(), copy); + return res; + } +} + +/** + * add all faces of top_face in simplices_per_dimension + */ +template<typename SimplexHandle> +void register_faces(std::vector< std::set<SimplexHandle> >& simplices_per_dimension, + const SimplexHandle& top_face) { + std::list<SimplexHandle> subfaces_list = subfaces(top_face); + for (auto& simplex : subfaces_list) { + simplices_per_dimension[simplex.dimension()].insert(simplex); + } +} + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_INTERNAL_TOP_FACES_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Trie.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Trie.h new file mode 100644 index 00000000..a43fa034 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/internal/Trie.h @@ -0,0 +1,256 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_INTERNAL_TRIE_H_ +#define SKELETON_BLOCKER_INTERNAL_TRIE_H_ + +#include <memory> +#include <vector> +#include <deque> +#include <set> + +namespace Gudhi { + +namespace skeleton_blocker { + +template<typename SimplexHandle> +struct Trie { + typedef SimplexHandle Simplex; + typedef typename SimplexHandle::Vertex_handle Vertex_handle; + + Vertex_handle v; + std::vector<std::shared_ptr<Trie> > childs; + // std::vector<std::unique_ptr<Trie> > childs; -> use of deleted function + private: + const Trie* parent_; + + public: + Trie() : parent_(0) { } + + Trie(Vertex_handle v_) : v(v_), parent_(0) { } + + Trie(Vertex_handle v_, Trie* parent) : v(v_), parent_(parent) { } + + bool operator==(const Trie& other) const { + return (v == other.v); + } + + void add_child(Trie* child) { + if (child) { + std::shared_ptr<Trie> ptr_to_add(child); + childs.push_back(ptr_to_add); + child->parent_ = this; + } + } + + typedef typename Simplex::Simplex_vertex_const_iterator Simplex_vertex_const_iterator; + + Trie* make_trie(Simplex_vertex_const_iterator s_it, Simplex_vertex_const_iterator s_end) { + if (s_it == s_end) { + return 0; + } else { + Trie* res = new Trie(*s_it); + Trie* child = make_trie(++s_it, s_end); + res->add_child(child); + return res; + } + } + + private: + // go down recursively in the tree while advancing the simplex iterator. + // when it reaches a leaf, it inserts the remaining that is not present + void add_simplex_helper(Simplex_vertex_const_iterator s_it, Simplex_vertex_const_iterator s_end) { + assert(*s_it == v); + ++s_it; + if (s_it == s_end) return; + if (!is_leaf()) { + for (auto child : childs) { + if (child->v == *s_it) + return child->add_simplex_helper(s_it, s_end); + } + // s_it is not found and needs to be inserted + } + // not leaf -> remaining of s needs to be inserted + Trie * son_with_what_remains_of_s(make_trie(s_it, s_end)); + add_child(son_with_what_remains_of_s); + return; + } + + void maximal_faces_helper(std::vector<Simplex>& res) const { + if (is_leaf()) res.push_back(simplex()); + else + for (auto child : childs) + child->maximal_faces_helper(res); + } + + public: + /** + * adds the simplex to the trie + */ + void add_simplex(const Simplex& s) { + if (s.empty()) return; + assert(v == s.first_vertex()); + add_simplex_helper(s.begin(), s.end()); + } + + std::vector<Simplex> maximal_faces() const { + std::vector<Simplex> res; + maximal_faces_helper(res); + return res; + } + + /** + * Goes to the root in the trie to consitute simplex + */ + void add_vertices_up_to_the_root(Simplex& res) const { + res.add_vertex(v); + if (parent_) + parent_->add_vertices_up_to_the_root(res); + } + + Simplex simplex() const { + Simplex res; + add_vertices_up_to_the_root(res); + return res; + } + + bool is_leaf() const { + return childs.empty(); + } + + bool is_root() const { + return parent_ == 0; + } + + const Trie* parent() { + return parent_; + } + + void remove_leaf() { + assert(is_leaf()); + if (!is_root()) + parent_->childs.erase(this); + } + + /** + * true iff the simplex corresponds to one node in the trie + */ + bool contains(const Simplex& s) const { + Trie const* current = this; + if (s.empty()) return true; + if (current->v != s.first_vertex()) return false; + auto s_pos = s.begin(); + ++s_pos; + while (s_pos != s.end() && current != 0) { + bool found = false; + for (const auto child : current->childs) { + if (child->v == *s_pos) { + ++s_pos; + current = child.get(); + found = true; + break; + } + } + if (!found) return false; + } + return current != 0; + } + + Trie* go_bottom_left() { + if (is_leaf()) + return this; + else + return (*childs.begin())->go_bottom_left(); + } + + friend std::ostream& operator<<(std::ostream& stream, const Trie& trie) { + stream << "T( " << trie.v << " "; + for (auto t : trie.childs) + stream << *t; + stream << ")"; + return stream; + } +}; + +template<typename SimplexHandle> +struct Tries { + typedef typename SimplexHandle::Vertex_handle Vertex_handle; + typedef SimplexHandle Simplex; + + typedef Trie<Simplex> STrie; + + template<typename SimpleHandleOutputIterator> + Tries(unsigned num_vertices, SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end) : + cofaces_(num_vertices, 0) { + for (auto i = 0u; i < num_vertices; ++i) + cofaces_[i] = new STrie(Vertex_handle(i)); + for (auto s_it = simplex_begin; s_it != simplex_end; ++s_it) { + if (s_it->dimension() >= 1) + cofaces_[s_it->first_vertex()]->add_simplex(*s_it); + } + } + + ~Tries() { + for (STrie* t : cofaces_) + delete t; + } + + // return a simplex that consists in all u such uv is an edge and u>v + + Simplex positive_neighbors(Vertex_handle v) const { + Simplex res; + for (auto child : cofaces_[v]->childs) + res.add_vertex(child->v); + return res; + } + + bool contains(const Simplex& s) const { + auto first_v = s.first_vertex(); + return cofaces_[first_v]->contains(s); + } + + friend std::ostream& operator<<(std::ostream& stream, const Tries& tries) { + for (auto trie : tries.cofaces_) + stream << *trie << std::endl; + return stream; + } + + // init_next_dimension must be called first + + std::vector<Simplex> next_dimension_simplices() const { + std::vector<Simplex> res; + while (!(to_see_.empty()) && (to_see_.front()->simplex().dimension() == current_dimension_)) { + res.emplace_back(to_see_.front()->simplex()); + for (auto child : to_see_.front()->childs) + to_see_.push_back(child.get()); + to_see_.pop_front(); + } + ++current_dimension_; + return res; + } + + void init_next_dimension() const { + for (auto trie : cofaces_) + to_see_.push_back(trie); + } + + private: + mutable std::deque<STrie*> to_see_; + mutable int current_dimension_ = 0; + std::vector<STrie*> cofaces_; +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_INTERNAL_TRIE_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h new file mode 100644 index 00000000..4f51f572 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h @@ -0,0 +1,122 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_BLOCKERS_ITERATORS_H_ +#define SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_BLOCKERS_ITERATORS_H_ + +#include <boost/iterator/iterator_facade.hpp> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + * @brief Iterator through the blockers of a vertex. + */ +// ReturnType = const Simplex* or Simplex* +// MapIteratorType = BlockerMapConstIterator or BlockerMapIterator + +template<typename MapIteratorType, typename ReturnType> +class Blocker_iterator_internal : public boost::iterator_facade< +Blocker_iterator_internal<MapIteratorType, ReturnType>, +ReturnType, +boost::forward_traversal_tag, +ReturnType +> { + private: + MapIteratorType current_position; + MapIteratorType end_of_map; + + public: + Blocker_iterator_internal() : current_position() { } + + Blocker_iterator_internal(MapIteratorType position, MapIteratorType end_of_map_) : + current_position(position), end_of_map(end_of_map_) { } + + bool equal(const Blocker_iterator_internal& other) const { + return current_position == other.current_position; + } + + void increment() { + goto_next_blocker(); + } + + ReturnType dereference() const { + return (current_position->second); + } + + private: + /** + * Let the current pair be (v,sigma) where v is a vertex and sigma is a blocker. + * If v is not the first vertex of sigma then we already have seen sigma as a blocker + * and we look for the next one. + */ + void goto_next_blocker() { + do { + ++current_position; + } while (!(current_position == end_of_map) && !first_time_blocker_is_seen()); + } + + bool first_time_blocker_is_seen() const { + return current_position->first == current_position->second->first_vertex(); + } +}; + +/** + * @brief Iterator through the blockers of a vertex + */ +// ReturnType = const Simplex* or Simplex* +// MapIteratorType = BlockerMapConstIterator or BlockerMapIterator + +template<typename MapIteratorType, typename ReturnType> +class Blocker_iterator_around_vertex_internal : public boost::iterator_facade< +Blocker_iterator_around_vertex_internal<MapIteratorType, ReturnType>, +ReturnType, +boost::forward_traversal_tag, +ReturnType +> { + private: + MapIteratorType current_position_; + + public: + Blocker_iterator_around_vertex_internal() : current_position_() { } + + Blocker_iterator_around_vertex_internal(MapIteratorType position) : + current_position_(position) { } + + Blocker_iterator_around_vertex_internal& operator=(Blocker_iterator_around_vertex_internal other) { + this->current_position_ = other.current_position_; + return *this; + } + + bool equal(const Blocker_iterator_around_vertex_internal& other) const { + return current_position_ == other.current_position_; + } + + void increment() { + current_position_++; + } + + ReturnType dereference() const { + return (current_position_->second); + } + + MapIteratorType current_position() { + return this->current_position_; + } +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_BLOCKERS_ITERATORS_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h new file mode 100644 index 00000000..154388a1 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h @@ -0,0 +1,135 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_EDGES_ITERATORS_H_ +#define SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_EDGES_ITERATORS_H_ + +#include <boost/iterator/iterator_facade.hpp> +#include <boost/graph/adjacency_list.hpp> + +#include <utility> // for pair<> + +namespace Gudhi { + +namespace skeleton_blocker { + +template<typename SkeletonBlockerComplex> +class Edge_around_vertex_iterator : public boost::iterator_facade <Edge_around_vertex_iterator<SkeletonBlockerComplex> +, typename SkeletonBlockerComplex::Edge_handle const, boost::forward_traversal_tag +, typename SkeletonBlockerComplex::Edge_handle const> { + friend class boost::iterator_core_access; + + typedef SkeletonBlockerComplex Complex; + typedef typename Complex::boost_adjacency_iterator boost_adjacency_iterator; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Edge_handle Edge_handle; + + private: + const Complex* complex; + Vertex_handle v; + + boost_adjacency_iterator current_; + boost_adjacency_iterator end_; + + public: + Edge_around_vertex_iterator() : complex(NULL) { } + + Edge_around_vertex_iterator(const Complex* complex_, Vertex_handle v_) : + complex(complex_), + v(v_) { + tie(current_, end_) = adjacent_vertices(v.vertex, complex->skeleton); + } + + /** + * returns an iterator to the end + */ + Edge_around_vertex_iterator(const Complex* complex_, Vertex_handle v_, int end) : + complex(complex_), + v(v_) { + tie(current_, end_) = adjacent_vertices(v.vertex, complex->skeleton); + set_end(); + } + + bool equal(const Edge_around_vertex_iterator& other) const { + return (complex == other.complex) && (v == other.v) && (current_ == other.current_) && (end_ == other.end_); + } + + void increment() { + if (current_ != end_) + ++current_; + } + + Edge_handle dereference() const { + return *(*complex)[std::make_pair(v, static_cast<Vertex_handle> (*current_))]; + } + + private: + // remove this ugly hack + void set_end() { + current_ = end_; + } +}; + +/** + *@brief Iterator on the edges of a simplicial complex. + * + */ +template<typename SkeletonBlockerComplex> +class Edge_iterator : public boost::iterator_facade <Edge_iterator<SkeletonBlockerComplex> +, typename SkeletonBlockerComplex::Edge_handle const +, boost::forward_traversal_tag +, typename SkeletonBlockerComplex::Edge_handle const> { + friend class boost::iterator_core_access; + + public: + typedef SkeletonBlockerComplex Complex; + typedef typename Complex::boost_edge_iterator boost_edge_iterator; + typedef typename Complex::Edge_handle Edge_handle; + + const Complex* complex; + std::pair<boost_edge_iterator, boost_edge_iterator> edge_iterator; + + Edge_iterator() : complex(NULL) { } + + Edge_iterator(const SkeletonBlockerComplex* complex_) : + complex(complex_), + edge_iterator(boost::edges(complex_->skeleton)) { } + + /** + * return an iterator to the end + */ + Edge_iterator(const SkeletonBlockerComplex* complex_, int end) : + complex(complex_), + edge_iterator(boost::edges(complex_->skeleton)) { + edge_iterator.first = edge_iterator.second; + } + + bool equal(const Edge_iterator& other) const { + return (complex == other.complex) && (edge_iterator == other.edge_iterator); + } + + void increment() { + if (edge_iterator.first != edge_iterator.second) { + ++(edge_iterator.first); + } + } + + Edge_handle dereference() const { + return (*(edge_iterator.first)); + } +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_EDGES_ITERATORS_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h new file mode 100644 index 00000000..7b43e05f --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h @@ -0,0 +1,20 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_ITERATORS_H_ +#define SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_ITERATORS_H_ + +#include <gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h> +#include <gudhi/Skeleton_blocker/iterators/Skeleton_blockers_edges_iterators.h> +#include <gudhi/Skeleton_blocker/iterators/Skeleton_blockers_blockers_iterators.h> +#include <gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h> +#include <gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h> + +#endif // SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_ITERATORS_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h new file mode 100644 index 00000000..920f8cb6 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_simplices_iterators.h @@ -0,0 +1,390 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_SIMPLICES_ITERATORS_H_ +#define SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_SIMPLICES_ITERATORS_H_ + +#include <gudhi/Skeleton_blocker_link_complex.h> +#include <gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h> +#include <gudhi/Skeleton_blocker/internal/Trie.h> +#include <gudhi/Debug_utils.h> + +#include <boost/iterator/iterator_facade.hpp> + +#include <memory> +#include <list> +#include <iostream> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + * Link may be Skeleton_blocker_link_complex<SkeletonBlockerComplex> to iterate over all + * simplices around a vertex OR + * Skeleton_blocker_superior_link_complex<SkeletonBlockerComplex> to iterate over all + * superior vertices around a vertex. + * The iteration is done by computing a trie with the link and doing a breadth-first traversal + * of the trie. + */ +template<typename SkeletonBlockerComplex, typename Link> +class Simplex_around_vertex_iterator : +public boost::iterator_facade < Simplex_around_vertex_iterator<SkeletonBlockerComplex, Link> +, typename SkeletonBlockerComplex::Simplex +, boost::forward_traversal_tag +, typename SkeletonBlockerComplex::Simplex +> { + friend class boost::iterator_core_access; + typedef SkeletonBlockerComplex Complex; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Edge_handle Edge_handle; + typedef typename Complex::Simplex Simplex; + + // Link_vertex_handle == Complex_Vertex_handle but this renaming helps avoiding confusion + typedef typename Link::Vertex_handle Link_vertex_handle; + + typedef typename Gudhi::skeleton_blocker::Trie<Simplex> Trie; + + private: + const Complex* complex; + Vertex_handle v; + std::shared_ptr<Link> link_v; + std::shared_ptr<Trie> trie; + // TODO(DS): use a deque instead + std::list<Trie*> nodes_to_be_seen; + + public: + Simplex_around_vertex_iterator() : complex(0) { } + + Simplex_around_vertex_iterator(const Complex* complex_, Vertex_handle v_) : + complex(complex_), + v(v_), + link_v(new Link(*complex_, v_)), + trie(new Trie(v_)) { + compute_trie_and_nodes_to_be_seen(); + } + + // TODO(DS): avoid useless copy + // TODO(DS): currently just work if copy begin iterator + Simplex_around_vertex_iterator(const Simplex_around_vertex_iterator& other) : + complex(other.complex), + v(other.v), + link_v(other.link_v), + trie(other.trie), + nodes_to_be_seen(other.nodes_to_be_seen) { + if (!other.is_end()) { + } + } + + /** + * returns an iterator to the end + */ + Simplex_around_vertex_iterator(const Complex* complex_, Vertex_handle v_, bool end) : + complex(complex_), + v(v_) { + set_end(); + } + + private: + void compute_trie_and_nodes_to_be_seen() { + // once we go through every simplices passing through v0 + // we remove v0. That way, it prevents from passing a lot of times + // though edges leaving v0. + // another solution would have been to provides an adjacency iterator + // to superior vertices that avoids lower ones. + while (!link_v->empty()) { + auto v0 = *(link_v->vertex_range().begin()); + trie->add_child(build_trie(v0, trie.get())); + link_v->remove_vertex(v0); + } + nodes_to_be_seen.push_back(trie.get()); + } + + Trie* build_trie(Link_vertex_handle link_vh, Trie* parent) { + Trie* res = new Trie(parent_vertex(link_vh), parent); + for (Link_vertex_handle nv : link_v->vertex_range(link_vh)) { + if (link_vh < nv) { + Simplex simplex_node_plus_nv(res->simplex()); + simplex_node_plus_nv.add_vertex(parent_vertex(nv)); + if (complex->contains(simplex_node_plus_nv)) { + res->add_child(build_trie(nv, res)); + } + } + } + return res; + } + + bool is_node_in_complex(Trie* trie) { + return true; + } + + Vertex_handle parent_vertex(Link_vertex_handle link_vh) const { + return complex->convert_handle_from_another_complex(*link_v, link_vh); + } + + public: + friend std::ostream& operator<<(std::ostream& stream, const Simplex_around_vertex_iterator& savi) { + stream << savi.trie << std::endl; + stream << "(" << savi.nodes_to_be_seen.size() << ") nodes to see\n"; + return stream; + } + + bool equal(const Simplex_around_vertex_iterator& other) const { + bool same_complex = (complex == other.complex); + if (!same_complex) + return false; + + if (!(v == other.v)) + return false; + + bool both_empty = nodes_to_be_seen.empty() && other.nodes_to_be_seen.empty(); + if (both_empty) + return true; + + bool both_non_empty = !nodes_to_be_seen.empty() && !other.nodes_to_be_seen.empty(); + + // one is empty the other is not + if (!both_non_empty) return false; + + bool same_node = (**(nodes_to_be_seen.begin()) == **(other.nodes_to_be_seen.begin())); + return same_node; + } + + void increment() { + assert(!is_end()); + Trie* first_node = nodes_to_be_seen.front(); + + nodes_to_be_seen.pop_front(); + + for (auto childs : first_node->childs) { + nodes_to_be_seen.push_back(childs.get()); + } + } + + Simplex dereference() const { + assert(!nodes_to_be_seen.empty()); + Trie* first_node = nodes_to_be_seen.front(); + return first_node->simplex(); + } + + Simplex get_trie_address() const { + assert(!nodes_to_be_seen.empty()); + return nodes_to_be_seen.front(); + } + + private: + void set_end() { + nodes_to_be_seen.clear(); + } + + bool is_end() const { + return nodes_to_be_seen.empty(); + } +}; + +template<typename SkeletonBlockerComplex> +class Simplex_iterator : +public boost::iterator_facade < Simplex_iterator<SkeletonBlockerComplex> +, typename SkeletonBlockerComplex::Simplex +, boost::forward_traversal_tag +, typename SkeletonBlockerComplex::Simplex +> { + typedef Skeleton_blocker_link_superior<SkeletonBlockerComplex> Link; + + friend class boost::iterator_core_access; + + template<class SkBlDS> friend class Skeleton_blocker_complex; + + typedef SkeletonBlockerComplex Complex; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Edge_handle Edge_handle; + typedef typename Complex::Simplex Simplex; + typedef typename Complex::Complex_vertex_iterator Complex_vertex_iterator; + typedef typename Link::Vertex_handle Link_vertex_handle; + + private: + const Complex* complex_; + Complex_vertex_iterator current_vertex_; + + typedef Simplex_around_vertex_iterator<SkeletonBlockerComplex, Link> SAVI; + SAVI current_simplex_around_current_vertex_; + SAVI simplices_around_current_vertex_end_; + + public: + Simplex_iterator() : complex_(0) { } + + Simplex_iterator(const Complex* complex) : + complex_(complex), + current_vertex_(complex->vertex_range().begin()), + current_simplex_around_current_vertex_(complex, *current_vertex_), + simplices_around_current_vertex_end_(complex, *current_vertex_, true) { + // should not be called if the complex is empty + assert(!complex->empty()); + } + + private: + Simplex_iterator(const Complex* complex, bool end) : + complex_(complex) { + set_end(); + } + + public: + Simplex_iterator(const Simplex_iterator& other) + : + complex_(other.complex_), + current_vertex_(other.current_vertex_), + current_simplex_around_current_vertex_(other.current_simplex_around_current_vertex_), + simplices_around_current_vertex_end_(other.simplices_around_current_vertex_end_) { } + + friend Simplex_iterator make_begin_iterator(const Complex* complex) { + if (complex->empty()) + return make_end_simplex_iterator(complex); + else + return Simplex_iterator(complex); + } + + friend Simplex_iterator make_end_simplex_iterator(const Complex* complex) { + return Simplex_iterator(complex, true); + } + + bool equal(const Simplex_iterator& other) const { + if (complex_ != other.complex_) return false; + if (current_vertex_ != other.current_vertex_) return false; + if (is_end() && other.is_end()) return true; + if (current_simplex_around_current_vertex_ != other.current_simplex_around_current_vertex_) + return false; + return true; + } + + void increment() { + if (current_simplex_around_current_vertex_ != simplices_around_current_vertex_end_) { + current_simplex_around_current_vertex_.increment(); + if (current_simplex_around_current_vertex_ == simplices_around_current_vertex_end_) + goto_next_vertex(); + } else { + goto_next_vertex(); + } + } + + void goto_next_vertex() { + current_vertex_.increment(); + if (!is_end()) { + current_simplex_around_current_vertex_ = SAVI(complex_, *current_vertex_); + simplices_around_current_vertex_end_ = SAVI(complex_, *current_vertex_, true); + } + } + + Simplex dereference() const { + return current_simplex_around_current_vertex_.dereference(); + } + + private: + void set_end() { + current_vertex_ = complex_->vertex_range().end(); + } + + bool is_end() const { + return (current_vertex_ == complex_->vertex_range().end()); + } +}; + +/** + * Iterator through the maximal faces of the coboundary of a simplex. + */ +template<typename SkeletonBlockerComplex, typename Link> +class Simplex_coboundary_iterator : +public boost::iterator_facade < Simplex_coboundary_iterator<SkeletonBlockerComplex, Link> +, typename SkeletonBlockerComplex::Simplex, boost::forward_traversal_tag, typename SkeletonBlockerComplex::Simplex> { + friend class boost::iterator_core_access; + typedef SkeletonBlockerComplex Complex; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Edge_handle Edge_handle; + typedef typename Complex::Simplex Simplex; + typedef typename Complex::Complex_vertex_iterator Complex_vertex_iterator; + + // Link_vertex_handle == Complex_Vertex_handle but this renaming helps avoiding confusion + typedef typename Link::Vertex_handle Link_vertex_handle; + + private: + const Complex* complex; + const Simplex& sigma; + std::shared_ptr<Link> link; + Complex_vertex_iterator current_vertex; + Complex_vertex_iterator link_vertex_end; + + public: + Simplex_coboundary_iterator() : complex(0) { } + + Simplex_coboundary_iterator(const Complex* complex_, const Simplex& sigma_) : + complex(complex_), + sigma(sigma_), + // need only vertices of the link hence the true flag + link(new Link(*complex_, sigma_, false, true)) { + auto link_vertex_range = link->vertex_range(); + current_vertex = link_vertex_range.begin(); + link_vertex_end = link_vertex_range.end(); + } + + Simplex_coboundary_iterator(const Simplex_coboundary_iterator& other) : + complex(other.complex), + sigma(other.sigma), + link(other.link), + current_vertex(other.current_vertex), + link_vertex_end(other.link_vertex_end) { } + + // returns an iterator to the end + Simplex_coboundary_iterator(const Complex* complex_, const Simplex& sigma_, bool end) : + complex(complex_), + sigma(sigma_) { + // to represent an end iterator without having to build a useless link, we use + // the convection that link is not initialized. + } + + private: + Vertex_handle parent_vertex(Link_vertex_handle link_vh) const { + return complex->convert_handle_from_another_complex(*link, link_vh); + } + + public: + friend std::ostream& operator<<(std::ostream& stream, const Simplex_coboundary_iterator& sci) { + return stream; + } + + // assume that iterator points to the same complex and comes from the same simplex + bool equal(const Simplex_coboundary_iterator& other) const { + assert(complex == other.complex && sigma == other.sigma); + if (is_end()) return other.is_end(); + if (other.is_end()) return is_end(); + return *current_vertex == *(other.current_vertex); + } + + void increment() { + ++current_vertex; + } + + Simplex dereference() const { + Simplex res(sigma); + res.add_vertex(parent_vertex(*current_vertex)); + return res; + } + + private: + bool is_end() const { + return !link || current_vertex == link_vertex_end; + } +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_SIMPLICES_ITERATORS_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h new file mode 100644 index 00000000..37c0b4d3 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_triangles_iterators.h @@ -0,0 +1,214 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_TRIANGLES_ITERATORS_H_ +#define SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_TRIANGLES_ITERATORS_H_ + +#include <boost/iterator/iterator_facade.hpp> +#include <memory> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + * \brief Iterator over the triangles that are + * adjacent to a vertex of the simplicial complex. + * \remark Will be removed soon -> dont look + */ +template<typename Complex, typename LinkType> +class Triangle_around_vertex_iterator : public boost::iterator_facade +< Triangle_around_vertex_iterator <Complex, LinkType> +, typename Complex::Simplex const +, boost::forward_traversal_tag +, typename Complex::Simplex const> { + friend class boost::iterator_core_access; + template<typename T> friend class Triangle_iterator; + private: + typedef typename LinkType::Vertex_handle Vertex_handle; + typedef typename LinkType::Root_vertex_handle Root_vertex_handle; + typedef typename LinkType::Simplex Simplex; + typedef typename Complex::Complex_edge_iterator Complex_edge_iterator_; + + const Complex* complex_; + Vertex_handle v_; + std::shared_ptr<LinkType> link_; + Complex_edge_iterator_ current_edge_; + bool is_end_; + + public: + Triangle_around_vertex_iterator(const Complex* complex, Vertex_handle v) : + complex_(complex), v_(v), link_(new LinkType(*complex, v_)), + current_edge_(link_->edge_range().begin()), + is_end_(current_edge_ == link_->edge_range().end()) { } + + /** + * @brief ugly hack to get an iterator to the end + */ + Triangle_around_vertex_iterator(const Complex* complex, Vertex_handle v, bool is_end) : + complex_(complex), v_(v), link_(0), is_end_(true) { } + + /** + * @brief ugly hack to get an iterator to the end + */ + Triangle_around_vertex_iterator() : + complex_(0), v_(-1), link_(0), is_end_(true) { } + + Triangle_around_vertex_iterator(const Triangle_around_vertex_iterator& other) { + v_ = other.v_; + complex_ = other.complex_; + is_end_ = other.is_end_; + + if (!is_end_) { + link_ = other.link_; + current_edge_ = other.current_edge_; + } + } + + bool equal(const Triangle_around_vertex_iterator& other) const { + return (complex_ == other.complex_) && ((finished() && other.finished()) || current_edge_ == other.current_edge_); + } + + Simplex dereference() const { + Root_vertex_handle v1 = (*link_)[*current_edge_].first(); + Root_vertex_handle v2 = (*link_)[*current_edge_].second(); + return Simplex(v_, *(complex_->get_address(v1)), *(complex_->get_address(v2))); + } + + void increment() { + ++current_edge_; + } + + private: + bool finished() const { + return is_end_ || (current_edge_ == link_->edge_range().end()); + } +}; + +/** + * \brief Iterator over the triangles of the + * simplicial complex. + * \remark Will be removed soon -> dont look + * + */ +template<typename SkeletonBlockerComplex> +class Triangle_iterator : public boost::iterator_facade< +Triangle_iterator <SkeletonBlockerComplex>, +typename SkeletonBlockerComplex::Simplex const +, boost::forward_traversal_tag +, typename SkeletonBlockerComplex::Simplex const> { + friend class boost::iterator_core_access; + private: + typedef typename SkeletonBlockerComplex::Vertex_handle Vertex_handle; + typedef typename SkeletonBlockerComplex::Root_vertex_handle Root_vertex_handle; + typedef typename SkeletonBlockerComplex::Simplex Simplex; + typedef typename SkeletonBlockerComplex::Superior_triangle_around_vertex_iterator STAVI; + typedef typename SkeletonBlockerComplex::Complex_vertex_iterator Complex_vertex_iterator; + + const SkeletonBlockerComplex* complex_; + Complex_vertex_iterator current_vertex_; + STAVI current_triangle_; + bool is_end_; + + public: + /* + * @remark assume that the complex is non-empty + */ + Triangle_iterator(const SkeletonBlockerComplex* complex) : + complex_(complex), + current_vertex_(complex->vertex_range().begin()), + current_triangle_(complex, *current_vertex_), // this line is problematic is the complex is empty + is_end_(false) { + assert(!complex->empty()); + gotoFirstTriangle(); + } + + private: + // goto to the first triangle or to the end if none + void gotoFirstTriangle() { + if (!is_finished() && current_triangle_.finished()) { + goto_next_vertex(); + } + } + + public: + /** + * @brief ugly hack to get an iterator to the end + * @remark assume that the complex is non-empty + */ + Triangle_iterator(const SkeletonBlockerComplex* complex, bool is_end) : + complex_(complex), + current_vertex_(complex->vertex_range().end()), + current_triangle_(), // xxx this line is problematic is the complex is empty + is_end_(true) { } + + Triangle_iterator& operator=(const Triangle_iterator & other) { + complex_ = other.complex_; + Complex_vertex_iterator current_vertex_; + STAVI current_triangle_; + return *this; + } + + bool equal(const Triangle_iterator& other) const { + bool both_are_finished = is_finished() && other.is_finished(); + bool both_arent_finished = !is_finished() && !other.is_finished(); + // if the two iterators are not finished, they must have the same state + return (complex_ == other.complex_) && (both_are_finished || ((both_arent_finished) && + current_vertex_ == other.current_vertex_ && + current_triangle_ == other.current_triangle_)); + } + + Simplex dereference() const { + return *current_triangle_; + } + + private: + // goto the next vertex that has a triangle pending or the + // end vertex iterator if none exists + void goto_next_vertex() { + // we must have consume all triangles passing through the vertex + assert(current_triangle_.finished()); + // we must not be done + assert(!is_finished()); + + ++current_vertex_; + + if (!is_finished()) { + current_triangle_ = STAVI(complex_, *current_vertex_); + if (current_triangle_.finished()) + goto_next_vertex(); + } + } + + public: + void increment() { + if (!current_triangle_.finished()) { + ++current_triangle_; // problem here + if (current_triangle_.finished()) + goto_next_vertex(); + } else { + assert(!is_finished()); + goto_next_vertex(); + } + } + + private: + bool is_finished() const { + return is_end_ || current_vertex_ == complex_->vertex_range().end(); + } +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_TRIANGLES_ITERATORS_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h new file mode 100644 index 00000000..49e94256 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker/iterators/Skeleton_blockers_vertices_iterators.h @@ -0,0 +1,158 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_VERTICES_ITERATORS_H_ +#define SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_VERTICES_ITERATORS_H_ + +#include <boost/iterator/iterator_facade.hpp> + +#include <utility> // for pair<> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + *@brief Iterator on the vertices of a simplicial complex + *@remark The increment operator go to the next active vertex. + *@remark Incrementation increases Vertex_handle. + */ +template<typename SkeletonBlockerComplex> +class Vertex_iterator : public boost::iterator_facade< Vertex_iterator <SkeletonBlockerComplex> +, typename SkeletonBlockerComplex::Vertex_handle const +, boost::forward_traversal_tag +, typename SkeletonBlockerComplex::Vertex_handle const> { + friend class boost::iterator_core_access; + + typedef typename SkeletonBlockerComplex::boost_vertex_iterator boost_vertex_iterator; + typedef typename SkeletonBlockerComplex::Vertex_handle Vertex_handle; + private: + const SkeletonBlockerComplex* complex; + std::pair<boost_vertex_iterator, boost_vertex_iterator> vertexIterator; + + + public: + Vertex_iterator() : complex(NULL) { } + + Vertex_iterator(const SkeletonBlockerComplex* complex_) : + complex(complex_), + vertexIterator(vertices(complex_->skeleton)) { + if (!finished() && !is_active()) { + goto_next_valid(); + } + } + + /** + * return an iterator to the end. + */ + Vertex_iterator(const SkeletonBlockerComplex* complex_, int end) : + complex(complex_), vertexIterator(vertices(complex_->skeleton)) { + vertexIterator.first = vertexIterator.second; + } + + public: + void increment() { + goto_next_valid(); + } + + Vertex_handle dereference() const { + return (Vertex_handle(*(vertexIterator.first))); + } + + bool equal(const Vertex_iterator& other) const { + return vertexIterator == other.vertexIterator && complex == other.complex; + } + + bool operator<(const Vertex_iterator& other) const { + return dereference() < other.dereference(); + } + + private: + bool finished() const { + return vertexIterator.first == vertexIterator.second; + } + + void goto_next_valid() { + ++vertexIterator.first; + if (!finished() && !is_active()) { + goto_next_valid(); + } + } + + bool is_active() const { + return ((*complex)[Vertex_handle(*vertexIterator.first)]).is_active(); + } +}; + +template<typename SkeletonBlockerComplex> +class Neighbors_vertices_iterator : public boost::iterator_facade < Neighbors_vertices_iterator<SkeletonBlockerComplex> +, typename SkeletonBlockerComplex::Vertex_handle const +, boost::forward_traversal_tag +, typename SkeletonBlockerComplex::Vertex_handle const> { + friend class boost::iterator_core_access; + + typedef SkeletonBlockerComplex Complex; + typedef typename Complex::boost_adjacency_iterator boost_adjacency_iterator; + typedef typename Complex::Vertex_handle Vertex_handle; + typedef typename Complex::Edge_handle Edge_handle; + + private: + const Complex* complex; + Vertex_handle v; + + boost_adjacency_iterator current_; + boost_adjacency_iterator end_; + + public: + Neighbors_vertices_iterator() : complex(NULL) { } + + Neighbors_vertices_iterator(const Complex* complex_, Vertex_handle v_) : + complex(complex_), + v(v_) { + tie(current_, end_) = adjacent_vertices(v.vertex, complex->skeleton); + } + + /** + * returns an iterator to the end + */ + Neighbors_vertices_iterator(const Complex* complex_, Vertex_handle v_, int end) : + complex(complex_), + v(v_) { + tie(current_, end_) = adjacent_vertices(v.vertex, complex->skeleton); + set_end(); + } + + void increment() { + if (current_ != end_) + ++current_; + } + + Vertex_handle dereference() const { + return (Vertex_handle(*current_)); + } + + bool equal(const Neighbors_vertices_iterator& other) const { + return (complex == other.complex) && (v == other.v) && (current_ == other.current_) && (end_ == other.end_); + } + + private: + // TODO(DS): remove this ugly hack + void set_end() { + current_ = end_; + } +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_ITERATORS_SKELETON_BLOCKERS_VERTICES_ITERATORS_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h new file mode 100644 index 00000000..125c6387 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_complex.h @@ -0,0 +1,1593 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_COMPLEX_H_ +#define SKELETON_BLOCKER_COMPLEX_H_ + +#include <gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h> +#include <gudhi/Skeleton_blocker_link_complex.h> +#include <gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h> +#include <gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h> +#include <gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h> +#include <gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h> +#include <gudhi/Skeleton_blocker/internal/Top_faces.h> +#include <gudhi/Skeleton_blocker/internal/Trie.h> +#include <gudhi/Debug_utils.h> + +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> +#include <boost/iterator/transform_iterator.hpp> +#include <boost/range/adaptor/map.hpp> + +#include <iostream> +#include <fstream> +#include <sstream> +#include <memory> +#include <map> +#include <list> +#include <set> +#include <vector> +#include <string> +#include <algorithm> +#include <utility> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + *@class Skeleton_blocker_complex + *@brief Abstract Simplicial Complex represented with a skeleton/blockers pair. + *@ingroup skbl + */ +template<class SkeletonBlockerDS> +class Skeleton_blocker_complex { + template<class ComplexType> friend class Vertex_iterator; + template<class ComplexType> friend class Neighbors_vertices_iterator; + template<class ComplexType> friend class Edge_iterator; + template<class ComplexType> friend class Edge_around_vertex_iterator; + + template<class ComplexType> friend class Skeleton_blocker_link_complex; + template<class ComplexType> friend class Skeleton_blocker_link_superior; + template<class ComplexType> friend class Skeleton_blocker_sub_complex; + + public: + /** + * @brief The type of stored vertex node, specified by the template SkeletonBlockerDS + */ + typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex; + + /** + * @brief The type of stored edge node, specified by the template SkeletonBlockerDS + */ + typedef typename SkeletonBlockerDS::Graph_edge Graph_edge; + + typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle; + + /** + * @brief The type of an handle to a vertex of the complex. + */ + typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle; + typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle; + + /** + * @brief A ordered set of integers that represents a simplex. + */ + typedef Skeleton_blocker_simplex<Vertex_handle> Simplex; + typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle; + + /** + * @brief Handle to a blocker of the complex. + */ + typedef Simplex* Blocker_handle; + + typedef typename Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_iterator; + typedef typename Simplex::Simplex_vertex_const_iterator Simplex_handle_iterator; + + protected: + typedef typename boost::adjacency_list<boost::setS, // edges + boost::vecS, // vertices + boost::undirectedS, Graph_vertex, Graph_edge> Graph; + // todo/remark : edges are not sorted, it heavily penalizes computation for SuperiorLink + // (eg Link with greater vertices) + // that burdens simplex iteration / complex initialization via list of simplices. + // to avoid that, one should modify the graph by storing two lists of adjacency for every + // vertex, the one with superior and the one with lower vertices, that way there is + // no more extra cost for computation of SuperiorLink + typedef typename boost::graph_traits<Graph>::vertex_iterator boost_vertex_iterator; + typedef typename boost::graph_traits<Graph>::edge_iterator boost_edge_iterator; + + protected: + typedef typename boost::graph_traits<Graph>::adjacency_iterator boost_adjacency_iterator; + + public: + /** + * @brief Handle to an edge of the complex. + */ + typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle; + + protected: + typedef std::multimap<Vertex_handle, Simplex *> BlockerMap; + typedef typename std::multimap<Vertex_handle, Simplex *>::value_type BlockerPair; + typedef typename std::multimap<Vertex_handle, Simplex *>::iterator BlockerMapIterator; + typedef typename std::multimap<Vertex_handle, Simplex *>::const_iterator BlockerMapConstIterator; + + protected: + size_t num_vertices_; + size_t num_blockers_; + + typedef Skeleton_blocker_complex_visitor<Vertex_handle> Visitor; + // typedef Visitor* Visitor_ptr; + Visitor* visitor; + + /** + * @details If 'x' is a Vertex_handle of a vertex in the complex then degree[x] = d is its degree. + * + * This quantity is updated when adding/removing edge. + * + * This is useful because the operation + * list.size() is done in linear time. + */ + std::vector<boost_vertex_handle> degree_; + Graph skeleton; /** 1-skeleton of the simplicial complex. */ + + /** Each vertex can access to the blockers passing through it. */ + BlockerMap blocker_map_; + + public: + ///////////////////////////////////////////////////////////////////////////// + /** @name Constructors, Destructors + */ + //@{ + + /** + *@brief constructs a simplicial complex with a given number of vertices and a visitor. + */ + explicit Skeleton_blocker_complex(size_t num_vertices_ = 0, Visitor* visitor_ = NULL) + : visitor(visitor_) { + clear(); + for (size_t i = 0; i < num_vertices_; ++i) { + add_vertex(); + } + } + + private: + // typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie; + typedef Trie<Simplex> STrie; + + public: + /** + * @brief Constructor with a list of simplices. + * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not). + */ + template<typename SimpleHandleOutputIterator> + Skeleton_blocker_complex(SimpleHandleOutputIterator simplices_begin, SimpleHandleOutputIterator simplices_end, + bool is_flag_complex = false, Visitor* visitor_ = NULL) + : num_vertices_(0), + num_blockers_(0), + visitor(visitor_) { + add_vertices_and_edges(simplices_begin, simplices_end); + + if (!is_flag_complex) + // need to compute blockers + add_blockers(simplices_begin, simplices_end); + } + + private: + /** + * Add vertices and edges of a simplex in one pass + */ + template<typename SimpleHandleOutputIterator> + void add_vertices_and_edges(SimpleHandleOutputIterator simplices_begin, SimpleHandleOutputIterator simplices_end) { + std::vector<std::pair<Vertex_handle, Vertex_handle>> edges; + // first pass to add vertices and edges + int num_vertex = -1; + for (auto s_it = simplices_begin; s_it != simplices_end; ++s_it) { + if (s_it->dimension() == 0) num_vertex = (std::max)(num_vertex, s_it->first_vertex().vertex); + if (s_it->dimension() == 1) edges.emplace_back(s_it->first_vertex(), s_it->last_vertex()); + } + while (num_vertex-- >= 0) add_vertex(); + + for (const auto& e : edges) + add_edge_without_blockers(e.first, e.second); + } + + template<typename SimpleHandleOutputIterator> + void add_blockers(SimpleHandleOutputIterator simplices_begin, SimpleHandleOutputIterator simplices_end) { + Tries<Simplex> tries(num_vertices(), simplices_begin, simplices_end); + tries.init_next_dimension(); + auto simplices(tries.next_dimension_simplices()); + + while (!simplices.empty()) { + simplices = tries.next_dimension_simplices(); + for (auto& sigma : simplices) { + // common_positive_neighbors is the set of vertices u such that + // for all s in sigma, us is an edge and u>s + Simplex common_positive_neighbors(tries.positive_neighbors(sigma.last_vertex())); + for (auto sigma_it = sigma.rbegin(); sigma_it != sigma.rend(); ++sigma_it) + if (sigma_it != sigma.rbegin()) + common_positive_neighbors.intersection(tries.positive_neighbors(*sigma_it)); + + for (auto x : common_positive_neighbors) { + // first test that all edges sx are here for all s in sigma + bool all_edges_here = true; + for (auto s : sigma) + if (!contains_edge(x, s)) { + all_edges_here = false; + break; + } + if (!all_edges_here) continue; + + // all edges of sigma \cup x are here + // we have a blocker if all proper faces of sigma \cup x + // are in the complex and if sigma \cup x is not in the complex + // the first is equivalent at checking if blocks(sigma \cup x) is true + // as all blockers of lower dimension have already been computed + sigma.add_vertex(x); + if (!tries.contains(sigma) && !blocks(sigma)) + add_blocker(sigma); + sigma.remove_vertex(x); + } + } + } + } + + public: + // We cannot use the default copy constructor since we need + // to make a copy of each of the blockers + + Skeleton_blocker_complex(const Skeleton_blocker_complex& copy) { + visitor = NULL; + degree_ = copy.degree_; + skeleton = Graph(copy.skeleton); + num_vertices_ = copy.num_vertices_; + + num_blockers_ = 0; + // we copy the blockers + for (auto blocker : copy.const_blocker_range()) { + add_blocker(*blocker); + } + } + + /** + */ + Skeleton_blocker_complex& operator=(const Skeleton_blocker_complex& copy) { + clear(); + visitor = NULL; + degree_ = copy.degree_; + skeleton = Graph(copy.skeleton); + num_vertices_ = copy.num_vertices_; + + num_blockers_ = 0; + // we copy the blockers + for (auto blocker : copy.const_blocker_range()) + add_blocker(*blocker); + return *this; + } + + /** + * return true if both complexes have the same simplices. + */ + bool operator==(const Skeleton_blocker_complex& other) const { + if (other.num_vertices() != num_vertices()) return false; + if (other.num_edges() != num_edges()) return false; + if (other.num_blockers() != num_blockers()) return false; + + for (auto v : vertex_range()) + if (!other.contains_vertex(v)) return false; + + for (auto e : edge_range()) + if (!other.contains_edge(first_vertex(e), second_vertex(e))) return false; + + for (const auto b : const_blocker_range()) + if (!other.contains_blocker(*b)) return false; + + return true; + } + + bool operator!=(const Skeleton_blocker_complex& other) const { + return !(*this == other); + } + + /** + * The destructor delete all blockers allocated. + */ + virtual ~Skeleton_blocker_complex() { + clear(); + } + + /** + * @details Clears the simplicial complex. After a call to this function, + * blockers are destroyed. The 1-skeleton and the set of blockers + * are both empty. + */ + virtual void clear() { + // xxx for now the responsabilty of freeing the visitor is for + // the user + visitor = NULL; + + degree_.clear(); + num_vertices_ = 0; + + remove_blockers(); + + skeleton.clear(); + } + + /** + *@brief allows to change the visitor. + */ + void set_visitor(Visitor* other_visitor) { + visitor = other_visitor; + } + + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Vertices operations + */ + //@{ + public: + /** + * @brief Return a local Vertex_handle of a vertex given a global one. + * @remark Assume that the vertex is present in the complex. + */ + Vertex_handle operator[](Root_vertex_handle global) const { + auto local(get_address(global)); + assert(local); + return *local; + } + + /** + * @brief Return the vertex node associated to local Vertex_handle. + * @remark Assume that the vertex is present in the complex. + */ + Graph_vertex& operator[](Vertex_handle address) { + assert( + 0 <= address.vertex && address.vertex < boost::num_vertices(skeleton)); + return skeleton[address.vertex]; + } + + /** + * @brief Return the vertex node associated to local Vertex_handle. + * @remark Assume that the vertex is present in the complex. + */ + const Graph_vertex& operator[](Vertex_handle address) const { + assert( + 0 <= address.vertex && address.vertex < boost::num_vertices(skeleton)); + return skeleton[address.vertex]; + } + + /** + * @brief Adds a vertex to the simplicial complex and returns its Vertex_handle. + * @remark Vertex representation is contiguous. + */ + Vertex_handle add_vertex() { + Vertex_handle address(boost::add_vertex(skeleton)); + num_vertices_++; + (*this)[address].activate(); + // safe since we now that we are in the root complex and the field 'address' and 'id' + // are identical for every vertices + (*this)[address].set_id(Root_vertex_handle(address.vertex)); + degree_.push_back(0); + if (visitor) + visitor->on_add_vertex(address); + return address; + } + + /** + * @brief Remove a vertex from the simplicial complex + * @remark It just deactivates the vertex with a boolean flag but does not + * remove it from vertices from complexity issues. + */ + void remove_vertex(Vertex_handle address) { + assert(contains_vertex(address)); + // We remove b + boost::clear_vertex(address.vertex, skeleton); + (*this)[address].deactivate(); + num_vertices_--; + degree_[address.vertex] = -1; + if (visitor) + visitor->on_remove_vertex(address); + } + + /** + */ + bool contains_vertex(Vertex_handle u) const { + Vertex_handle num_vertices(boost::num_vertices(skeleton)); + if (u.vertex < 0 || u.vertex >= num_vertices) + return false; + return (*this)[u].is_active(); + } + + /** + */ + bool contains_vertex(Root_vertex_handle u) const { + boost::optional<Vertex_handle> address = get_address(u); + return address && (*this)[*address].is_active(); + } + + /** + * @return true iff the simplicial complex contains all vertices + * of simplex sigma + */ + bool contains_vertices(const Simplex & sigma) const { + for (auto vertex : sigma) + if (!contains_vertex(vertex)) + return false; + return true; + } + + /** + * @brief Given an Id return the address of the vertex having this Id in the complex. + * @remark For a simplicial complex, the address is the id but it may not be the case for a SubComplex. + */ + virtual boost::optional<Vertex_handle> get_address(Root_vertex_handle id) const { + boost::optional<Vertex_handle> res; + int num_vertices = boost::num_vertices(skeleton); + if (id.vertex < num_vertices) + res = Vertex_handle(id.vertex); + return res; + } + + /** + * return the id of a vertex of adress local present in the graph + */ + Root_vertex_handle get_id(Vertex_handle local) const { + assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton)); + return (*this)[local].get_id(); + } + + /** + * @brief Convert an address of a vertex of a complex to the address in + * the current complex. + * @details + * If the current complex is a sub (or sup) complex of 'other', it converts + * the address of a vertex v expressed in 'other' to the address of the vertex + * v in the current one. + * @remark this methods uses Root_vertex_handle to identify the vertex and + * assumes the vertex is present in the current complex. + */ + Vertex_handle convert_handle_from_another_complex(const Skeleton_blocker_complex& other, + Vertex_handle vh_in_other) const { + auto vh_in_current_complex = get_address(other.get_id(vh_in_other)); + assert(vh_in_current_complex); + return *vh_in_current_complex; + } + + /** + * @brief return the graph degree of a vertex. + */ + int degree(Vertex_handle local) const { + assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton)); + return degree_[local.vertex]; + } + + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Edges operations + */ + //@{ + public: + /** + * @brief return an edge handle if the two vertices forms + * an edge in the complex + */ + boost::optional<Edge_handle> operator[]( + const std::pair<Vertex_handle, Vertex_handle>& ab) const { + boost::optional<Edge_handle> res; + std::pair<Edge_handle, bool> edge_pair( + boost::edge(ab.first.vertex, ab.second.vertex, skeleton)); + if (edge_pair.second) + res = edge_pair.first; + return res; + } + + /** + * @brief returns the stored node associated to an edge + */ + Graph_edge& operator[](Edge_handle edge_handle) { + return skeleton[edge_handle]; + } + + /** + * @brief returns the stored node associated to an edge + */ + const Graph_edge& operator[](Edge_handle edge_handle) const { + return skeleton[edge_handle]; + } + + /** + * @brief returns the first vertex of an edge + * @details it assumes that the edge is present in the complex + */ + Vertex_handle first_vertex(Edge_handle edge_handle) const { + return static_cast<Vertex_handle> (source(edge_handle, skeleton)); + } + + /** + * @brief returns the first vertex of an edge + * @details it assumes that the edge is present in the complex + */ + Vertex_handle second_vertex(Edge_handle edge_handle) const { + return static_cast<Vertex_handle> (target(edge_handle, skeleton)); + } + + /** + * @brief returns the simplex made with the two vertices of the edge + * @details it assumes that the edge is present in the complex + + */ + Simplex get_vertices(Edge_handle edge_handle) const { + auto edge((*this)[edge_handle]); + return Simplex((*this)[edge.first()], (*this)[edge.second()]); + } + + /** + * @brief Adds an edge between vertices a and b. + * @details For instance, the complex contains edge 01 and 12, then calling + * add_edge with vertex 0 and 2 will create a complex containing + * the edges 01, 12, 20 but not the triangle 012 (and hence this complex + * will contains a blocker 012). + */ + Edge_handle add_edge(Vertex_handle a, Vertex_handle b) { + // if the edge is already there we musnt go further + // as we may add blockers that should not be here + if (contains_edge(a, b)) + return *((*this)[std::make_pair(a, b)]); + auto res = add_edge_without_blockers(a, b); + add_blockers_after_simplex_insertion(Simplex(a, b)); + return res; + } + + /** + * @brief Adds all edges of s in the complex. + */ + void add_edge(const Simplex& s) { + for (auto i = s.begin(); i != s.end(); ++i) + for (auto j = i; ++j != s.end(); /**/) + add_edge(*i, *j); + } + + /** + * @brief Adds an edge between vertices a and b without blockers. + * @details For instance, the complex contains edge 01 and 12, then calling + * add_edge with vertex 0 and 2 will create a complex containing + * the triangle 012. + */ + Edge_handle add_edge_without_blockers(Vertex_handle a, Vertex_handle b) { + assert(contains_vertex(a) && contains_vertex(b)); + assert(a != b); + + auto edge_handle((*this)[std::make_pair(a, b)]); + if (!edge_handle) { + edge_handle = boost::add_edge(a.vertex, b.vertex, skeleton).first; + (*this)[*edge_handle].setId(get_id(a), get_id(b)); + degree_[a.vertex]++; + degree_[b.vertex]++; + if (visitor) + visitor->on_add_edge_without_blockers(a, b); + } + return *edge_handle; + } + + /** + * @brief Adds all edges of s in the complex without adding blockers. + */ + void add_edge_without_blockers(Simplex s) { + for (auto i = s.begin(); i != s.end(); ++i) { + for (auto j = i; ++j != s.end(); /**/) + add_edge_without_blockers(*i, *j); + } + } + + /** + * @brief Removes an edge from the simplicial complex and all its cofaces. + * @details returns the former Edge_handle representing the edge + */ + virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b) { + bool found; + Edge_handle edge; + tie(edge, found) = boost::edge(a.vertex, b.vertex, skeleton); + if (found) { + if (visitor) + visitor->on_remove_edge(a, b); + boost::remove_edge(a.vertex, b.vertex, skeleton); + degree_[a.vertex]--; + degree_[b.vertex]--; + } + return edge; + } + + /** + * @brief Removes edge and its cofaces from the simplicial complex. + */ + void remove_edge(Edge_handle edge) { + assert(contains_vertex(first_vertex(edge))); + assert(contains_vertex(second_vertex(edge))); + remove_edge(first_vertex(edge), second_vertex(edge)); + } + + /** + * @brief The complex is reduced to its set of vertices. + * All the edges and blockers are removed. + */ + void keep_only_vertices() { + remove_blockers(); + + for (auto u : vertex_range()) { + while (this->degree(u) > 0) { + Vertex_handle v(*(adjacent_vertices(u.vertex, this->skeleton).first)); + this->remove_edge(u, v); + } + } + } + + /** + * @return true iff the simplicial complex contains an edge between + * vertices a and b + */ + bool contains_edge(Vertex_handle a, Vertex_handle b) const { + // if (a.vertex<0 || b.vertex <0) return false; + return boost::edge(a.vertex, b.vertex, skeleton).second; + } + + /** + * @return true iff the simplicial complex contains all vertices + * and all edges of simplex sigma + */ + bool contains_edges(const Simplex & sigma) const { + for (auto i = sigma.begin(); i != sigma.end(); ++i) { + if (!contains_vertex(*i)) + return false; + for (auto j = i; ++j != sigma.end();) { + if (!contains_edge(*i, *j)) + return false; + } + } + return true; + } + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Blockers operations + */ + //@{ + + /** + * @brief Adds the simplex to the set of blockers and + * returns a Blocker_handle toward it if was not present before and 0 otherwise. + */ + Blocker_handle add_blocker(const Simplex& blocker) { + assert(blocker.dimension() > 1); + if (contains_blocker(blocker)) { + return 0; + } else { + if (visitor) + visitor->on_add_blocker(blocker); + Blocker_handle blocker_pt = new Simplex(blocker); + num_blockers_++; + auto vertex = blocker_pt->begin(); + while (vertex != blocker_pt->end()) { + blocker_map_.insert(BlockerPair(*vertex, blocker_pt)); + ++vertex; + } + return blocker_pt; + } + } + + protected: + /** + * @brief Adds the simplex to the set of blockers + */ + void add_blocker(Blocker_handle blocker) { + if (contains_blocker(*blocker)) { + // std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl; + return; + } else { + if (visitor) + visitor->on_add_blocker(*blocker); + num_blockers_++; + auto vertex = blocker->begin(); + while (vertex != blocker->end()) { + blocker_map_.insert(BlockerPair(*vertex, blocker)); + ++vertex; + } + } + } + + protected: + /** + * Removes sigma from the blocker map of vertex v + */ + void remove_blocker(const Blocker_handle sigma, Vertex_handle v) { + Complex_blocker_around_vertex_iterator blocker; + for (blocker = blocker_range(v).begin(); blocker != blocker_range(v).end(); + ++blocker) { + if (*blocker == sigma) + break; + } + if (*blocker != sigma) { + std::cerr + << "bug ((*blocker).second == sigma) ie try to remove a blocker not present\n"; + assert(false); + } else { + blocker_map_.erase(blocker.current_position()); + } + } + + public: + /** + * @brief Removes the simplex from the set of blockers. + * @remark sigma has to belongs to the set of blockers + */ + void remove_blocker(const Blocker_handle sigma) { + for (auto vertex : *sigma) + remove_blocker(sigma, vertex); + num_blockers_--; + } + + /** + * @brief Remove all blockers, in other words, it expand the simplicial + * complex to the smallest flag complex that contains it. + */ + void remove_blockers() { + // Desallocate the blockers + while (!blocker_map_.empty()) { + delete_blocker(blocker_map_.begin()->second); + } + num_blockers_ = 0; + blocker_map_.clear(); + } + + protected: + /** + * Removes the simplex sigma from the set of blockers. + * sigma has to belongs to the set of blockers + * + * @remark contrarily to delete_blockers does not call the destructor + */ + void remove_blocker(const Simplex& sigma) { + assert(contains_blocker(sigma)); + for (auto vertex : sigma) + remove_blocker(sigma, vertex); + num_blockers_--; + } + + public: + /** + * Removes the simplex s from the set of blockers + * and desallocate s. + */ + void delete_blocker(Blocker_handle sigma) { + if (visitor) + visitor->on_delete_blocker(sigma); + remove_blocker(sigma); + delete sigma; + } + + /** + * @return true iff s is a blocker of the simplicial complex + */ + bool contains_blocker(const Blocker_handle s) const { + if (s->dimension() < 2) + return false; + + Vertex_handle a = s->first_vertex(); + + for (const auto blocker : const_blocker_range(a)) { + if (s == *blocker) + return true; + } + return false; + } + + /** + * @return true iff s is a blocker of the simplicial complex + */ + bool contains_blocker(const Simplex & s) const { + if (s.dimension() < 2) + return false; + + Vertex_handle a = s.first_vertex(); + + for (auto blocker : const_blocker_range(a)) { + if (s == *blocker) + return true; + } + return false; + } + + private: + /** + * @return true iff a blocker of the simplicial complex + * is a face of sigma. + */ + bool blocks(const Simplex & sigma) const { + for (auto s : sigma) + for (auto blocker : const_blocker_range(s)) + if (sigma.contains(*blocker)) + return true; + return false; + } + + //@} + + protected: + /** + * @details Adds to simplex the neighbours of v e.g. \f$ n \leftarrow n \cup N(v) \f$. + * If keep_only_superior is true then only vertices that are greater than v are added. + */ + virtual void add_neighbours(Vertex_handle v, Simplex & n, + bool keep_only_superior = false) const { + boost_adjacency_iterator ai, ai_end; + for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end; + ++ai) { + Vertex_handle value(*ai); + if (keep_only_superior) { + if (value > v.vertex) { + n.add_vertex(value); + } + } else { + n.add_vertex(value); + } + } + } + + /** + * @details Add to simplex res all vertices which are + * neighbours of alpha: ie \f$ res \leftarrow res \cup N(alpha) \f$. + * + * If 'keep_only_superior' is true then only vertices that are greater than alpha are added. + * todo revoir + * + */ + virtual void add_neighbours(const Simplex &alpha, Simplex & res, + bool keep_only_superior = false) const { + res.clear(); + auto alpha_vertex = alpha.begin(); + add_neighbours(*alpha_vertex, res, keep_only_superior); + for (alpha_vertex = (alpha.begin())++; alpha_vertex != alpha.end(); + ++alpha_vertex) + keep_neighbours(*alpha_vertex, res, keep_only_superior); + } + + /** + * @details Remove from simplex n all vertices which are + * not neighbours of v e.g. \f$ res \leftarrow res \cap N(v) \f$. + * If 'keep_only_superior' is true then only vertices that are greater than v are keeped. + */ + virtual void keep_neighbours(Vertex_handle v, Simplex& res, + bool keep_only_superior = false) const { + Simplex nv; + add_neighbours(v, nv, keep_only_superior); + res.intersection(nv); + } + + /** + * @details Remove from simplex all vertices which are + * neighbours of v eg \f$ res \leftarrow res \setminus N(v) \f$. + * If 'keep_only_superior' is true then only vertices that are greater than v are added. + */ + virtual void remove_neighbours(Vertex_handle v, Simplex & res, + bool keep_only_superior = false) const { + Simplex nv; + add_neighbours(v, nv, keep_only_superior); + res.difference(nv); + } + + public: + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex> Link_complex; + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Link_complex link(Vertex_handle v) const { + return Link_complex(*this, Simplex(v)); + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Link_complex link(Edge_handle edge) const { + return Link_complex(*this, edge); + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Link_complex link(const Simplex& simplex) const { + return Link_complex(*this, simplex); + } + + /** + * @brief Compute the local vertices of 's' in the current complex + * If one of them is not present in the complex then the return value is uninitialized. + * + * + */ + // xxx rename get_address et place un using dans sub_complex + + boost::optional<Simplex> get_simplex_address( + const Root_simplex_handle& s) const { + boost::optional<Simplex> res; + + Simplex s_address; + // Root_simplex_const_iterator i; + for (auto i = s.begin(); i != s.end(); ++i) { + boost::optional<Vertex_handle> address = get_address(*i); + if (!address) + return res; + else + s_address.add_vertex(*address); + } + res = s_address; + return res; + } + + /** + * @brief returns a simplex with vertices which are the id of vertices of the + * argument. + */ + Root_simplex_handle get_id(const Simplex& local_simplex) const { + Root_simplex_handle global_simplex; + for (auto x = local_simplex.begin(); x != local_simplex.end(); ++x) { + global_simplex.add_vertex(get_id(*x)); + } + return global_simplex; + } + + /** + * @brief returns true iff the simplex s belongs to the simplicial + * complex. + */ + virtual bool contains(const Simplex & s) const { + if (s.dimension() == -1) { + return false; + } else if (s.dimension() == 0) { + return contains_vertex(s.first_vertex()); + } else { + return (contains_edges(s) && !blocks(s)); + } + } + + /* + * @brief returnrs true iff the complex is empty. + */ + bool empty() const { + return num_vertices() == 0; + } + + /* + * @brief returns the number of vertices in the complex. + */ + int num_vertices() const { + // remark boost::num_vertices(skeleton) counts deactivated vertices + return num_vertices_; + } + + /* + * @brief returns the number of edges in the complex. + * @details currently in O(n) + */ + // todo cache the value + + int num_edges() const { + return boost::num_edges(skeleton); + } + + int num_triangles() const { + auto triangles = triangle_range(); + return std::distance(triangles.begin(), triangles.end()); + } + + /* + * @brief returns the number of simplices of a given dimension in the complex. + */ + size_t num_simplices() const { + auto simplices = complex_simplex_range(); + return std::distance(simplices.begin(), simplices.end()); + } + + /* + * @brief returns the number of simplices of a given dimension in the complex. + */ + size_t num_simplices(int dimension) const { + // TODO(DS): iterator on k-simplices + size_t res = 0; + for (const auto& s : complex_simplex_range()) + if (s.dimension() == dimension) + ++res; + return res; + } + + /* + * @brief returns the number of blockers in the complex. + */ + size_t num_blockers() const { + return num_blockers_; + } + + /* + * @brief returns true iff the graph of the 1-skeleton of the complex is complete. + */ + bool complete() const { + return (num_vertices() * (num_vertices() - 1)) / 2 == num_edges(); + } + + /** + * @brief returns the number of connected components in the graph of the 1-skeleton. + */ + int num_connected_components() const { + int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices(); + std::vector<int> component(skeleton.vertex_set().size()); + return boost::connected_components(this->skeleton, &component[0]) + - num_vert_collapsed; + } + + /** + * @brief %Test if the complex is a cone. + * @details Runs in O(n) where n is the number of vertices. + */ + bool is_cone() const { + if (num_vertices() == 0) + return false; + if (num_vertices() == 1) + return true; + for (auto vi : vertex_range()) { + // xxx todo faire une methode bool is_in_blocker(Vertex_handle) + if (blocker_map_.find(vi) == blocker_map_.end()) { + // no blocker passes through the vertex, we just need to + // check if the current vertex is linked to all others vertices of the complex + if (degree_[vi.vertex] == num_vertices() - 1) + return true; + } + } + return false; + } + + //@} + /** @name Simplification operations + */ + //@{ + + /** + * Returns true iff the blocker 'sigma' is popable. + * To define popable, let us call 'L' the complex that + * consists in the current complex without the blocker 'sigma'. + * A blocker 'sigma' is then "popable" if the link of 'sigma' + * in L is reducible. + * + */ + bool is_popable_blocker(Blocker_handle sigma) const; + + /** + * Removes all the popable blockers of the complex and delete them. + * @returns the number of popable blockers deleted + */ + void remove_popable_blockers(); + + /** + * Removes all the popable blockers of the complex passing through v and delete them. + */ + void remove_popable_blockers(Vertex_handle v); + + /** + * @brief Removes all the popable blockers of the complex passing through v and delete them. + * Also remove popable blockers in the neighborhood if they became popable. + * + */ + void remove_all_popable_blockers(Vertex_handle v); + + /** + * Remove the star of the vertex 'v' + */ + void remove_star(Vertex_handle v); + + private: + /** + * after removing the star of a simplex, blockers sigma that contains this simplex must be removed. + * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added + * whenever the dimension of tau is at least 2. + */ + void update_blockers_after_remove_star_of_vertex_or_edge(const Simplex& simplex_to_be_removed); + + public: + /** + * Remove the star of the edge connecting vertices a and b. + * @returns the number of blocker that have been removed + */ + void remove_star(Vertex_handle a, Vertex_handle b); + + /** + * Remove the star of the edge 'e'. + */ + void remove_star(Edge_handle e); + + /** + * Remove the star of the simplex 'sigma' which needs to belong to the complex + */ + void remove_star(const Simplex& sigma); + + /** + * @brief add a simplex and all its faces. + * @details the simplex must have dimension greater than one (otherwise use add_vertex or add_edge_without_blockers). + */ + void add_simplex(const Simplex& sigma); + + private: + void add_blockers_after_simplex_insertion(Simplex s); + + /** + * remove all blockers that contains sigma + */ + void remove_blocker_containing_simplex(const Simplex& sigma); + + /** + * remove all blockers that contains sigma + */ + void remove_blocker_include_in_simplex(const Simplex& sigma); + + public: + enum simplifiable_status { + NOT_HOMOTOPY_EQ, MAYBE_HOMOTOPY_EQ, HOMOTOPY_EQ + }; + + simplifiable_status is_remove_star_homotopy_preserving(const Simplex& simplex) { + // todo write a virtual method 'link' in Skeleton_blocker_complex which will be overloaded by the current one of + // Skeleton_blocker_geometric_complex + // then call it there to build the link and return the value of link.is_contractible() + return MAYBE_HOMOTOPY_EQ; + } + + enum contractible_status { + NOT_CONTRACTIBLE, MAYBE_CONTRACTIBLE, CONTRACTIBLE + }; + + /** + * @brief %Test if the complex is reducible using a strategy defined in the class + * (by default it tests if the complex is a cone) + * @details Note that NO could be returned if some invariant ensures that the complex + * is not a point (for instance if the euler characteristic is different from 1). + * This function will surely have to return MAYBE in some case because the + * associated problem is undecidable but it in practice, it can often + * be solved with the help of geometry. + */ + virtual contractible_status is_contractible() const { + if (this->is_cone()) { + return CONTRACTIBLE; + } else { + return MAYBE_CONTRACTIBLE; + } + } + //@} + + /** @name Edge contraction operations + */ + //@{ + + /** + * @return If ignore_popable_blockers is true + * then the result is true iff the link condition at edge ab is satisfied + * or equivalently iff no blocker contains ab. + * If ignore_popable_blockers is false then the + * result is true iff all blocker containing ab are popable. + */ + bool link_condition(Vertex_handle a, Vertex_handle b, bool ignore_popable_blockers = false) const { + for (auto blocker : this->const_blocker_range(a)) + if (blocker->contains(b)) { + // false if ignore_popable_blockers is false + // otherwise the blocker has to be popable + return ignore_popable_blockers && is_popable_blocker(blocker); + } + return true; + } + + /** + * @return If ignore_popable_blockers is true + * then the result is true iff the link condition at edge ab is satisfied + * or equivalently iff no blocker contains ab. + * If ignore_popable_blockers is false then the + * result is true iff all blocker containing ab are popable. + */ + bool link_condition(Edge_handle e, bool ignore_popable_blockers = false) const { + const Graph_edge& edge = (*this)[e]; + assert(this->get_address(edge.first())); + assert(this->get_address(edge.second())); + Vertex_handle a(*this->get_address(edge.first())); + Vertex_handle b(*this->get_address(edge.second())); + return link_condition(a, b, ignore_popable_blockers); + } + + protected: + /** + * Compute simplices beta such that a.beta is an order 0 blocker + * that may be used to construct a new blocker after contracting ab. + * It requires that the link condition is satisfied. + */ + void tip_blockers(Vertex_handle a, Vertex_handle b, std::vector<Simplex> & buffer) const; + + private: + /** + * @brief "Replace" the edge 'bx' by the edge 'ax'. + * Assume that the edge 'bx' was present whereas 'ax' was not. + * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'. + * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added + * and just before edge 'bx' had been removed. That way, it can + * eventually access to information of 'ax'. + */ + void swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x); + + private: + /** + * @brief removes all blockers passing through the edge 'ab' + */ + void delete_blockers_around_vertex(Vertex_handle v); + + /** + * @brief removes all blockers passing through the edge 'ab' + */ + void delete_blockers_around_edge(Vertex_handle a, Vertex_handle b); + + public: + /** + * Contracts the edge. + * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, + * it removes first all blockers passing through 'ab' + */ + void contract_edge(Edge_handle edge) { + contract_edge(this->first_vertex(edge), this->second_vertex(edge)); + } + + /** + * Contracts the edge connecting vertices a and b. + * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, + * it removes first all blockers passing through 'ab' + */ + void contract_edge(Vertex_handle a, Vertex_handle b); + + private: + void get_blockers_to_be_added_after_contraction(Vertex_handle a, Vertex_handle b, + std::set<Simplex>& blockers_to_add); + /** + * delete all blockers that passes through a or b + */ + void delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b); + void update_edges_after_contraction(Vertex_handle a, Vertex_handle b); + void notify_changed_edges(Vertex_handle a); + //@} + + public: + ///////////////////////////////////////////////////////////////////////////// + /** @name Vertex iterators + */ + //@{ + typedef Vertex_iterator<Skeleton_blocker_complex> Complex_vertex_iterator; + + // + // Range over the vertices of the simplicial complex. + // Methods .begin() and .end() return a Complex_vertex_iterator. + // + typedef boost::iterator_range<Complex_vertex_iterator> Complex_vertex_range; + + /** + * @brief Returns a Complex_vertex_range over all vertices of the complex + */ + Complex_vertex_range vertex_range() const { + auto begin = Complex_vertex_iterator(this); + auto end = Complex_vertex_iterator(this, 0); + return Complex_vertex_range(begin, end); + } + + typedef Neighbors_vertices_iterator<Skeleton_blocker_complex> Complex_neighbors_vertices_iterator; + + + typedef boost::iterator_range<Complex_neighbors_vertices_iterator> Complex_neighbors_vertices_range; + + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v + */ + Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const { + auto begin = Complex_neighbors_vertices_iterator(this, v); + auto end = Complex_neighbors_vertices_iterator(this, v, 0); + return Complex_neighbors_vertices_range(begin, end); + } + + //@} + + /** @name Edge iterators + */ + //@{ + + typedef Edge_iterator<Skeleton_blocker_complex> Complex_edge_iterator; + + + typedef boost::iterator_range<Complex_edge_iterator> Complex_edge_range; + + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex + */ + Complex_edge_range edge_range() const { + auto begin = Complex_edge_iterator(this); + auto end = Complex_edge_iterator(this, 0); + return Complex_edge_range(begin, end); + } + + + typedef Edge_around_vertex_iterator<Skeleton_blocker_complex> Complex_edge_around_vertex_iterator; + + + typedef boost::iterator_range <Complex_edge_around_vertex_iterator> Complex_edge_around_vertex_range; + + /** + * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes + * through 'v' + */ + Complex_edge_around_vertex_range edge_range(Vertex_handle v) const { + auto begin = Complex_edge_around_vertex_iterator(this, v); + auto end = Complex_edge_around_vertex_iterator(this, v, 0); + return Complex_edge_around_vertex_range(begin, end); + } + + //@} + + /** @name Triangles iterators + */ + //@{ + private: + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link; + typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link; + + public: + typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex, Superior_link> + Superior_triangle_around_vertex_iterator; + typedef boost::iterator_range < Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link> > + Complex_triangle_around_vertex_range; + + /** + * @brief Range over triangles around a vertex of the simplicial complex. + * Methods .begin() and .end() return a Triangle_around_vertex_iterator. + * + */ + Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const { + auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v); + auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex, Link>(this, v, 0); + return Complex_triangle_around_vertex_range(begin, end); + } + + typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range; + typedef Triangle_iterator<Skeleton_blocker_complex> Complex_triangle_iterator; + + /** + * @brief Range over triangles of the simplicial complex. + * Methods .begin() and .end() return a Triangle_around_vertex_iterator. + * + */ + Complex_triangle_range triangle_range() const { + auto end = Triangle_iterator<Skeleton_blocker_complex>(this, 0); + if (empty()) { + return Complex_triangle_range(end, end); + } else { + auto begin = Triangle_iterator<Skeleton_blocker_complex>(this); + return Complex_triangle_range(begin, end); + } + } + + //@} + + /** @name Simplices iterators + */ + //@{ + typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex, Link> Complex_simplex_around_vertex_iterator; + + /** + * @brief Range over the simplices of the simplicial complex around a vertex. + * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator. + */ + typedef boost::iterator_range < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range; + + /** + * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex + */ + Complex_simplex_around_vertex_range star_simplex_range(Vertex_handle v) const { + assert(contains_vertex(v)); + return Complex_simplex_around_vertex_range( + Complex_simplex_around_vertex_iterator(this, v), + Complex_simplex_around_vertex_iterator(this, v, true)); + } + typedef Simplex_coboundary_iterator<Skeleton_blocker_complex, Link> Complex_simplex_coboundary_iterator; + + /** + * @brief Range over the simplices of the coboundary of a simplex. + * Methods .begin() and .end() return a Complex_simplex_coboundary_iterator. + */ + typedef boost::iterator_range < Complex_simplex_coboundary_iterator > Complex_coboundary_range; + + /** + * @brief Returns a Complex_simplex_coboundary_iterator over the simplices of the coboundary of a simplex. + */ + Complex_coboundary_range coboundary_range(const Simplex& s) const { + assert(contains(s)); + return Complex_coboundary_range(Complex_simplex_coboundary_iterator(this, s), + Complex_simplex_coboundary_iterator(this, s, true)); + } + + // typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator; + typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator; + + typedef boost::iterator_range < Complex_simplex_iterator > Complex_simplex_range; + + /** + * @brief Returns a Complex_simplex_range over all the simplices of the complex + */ + Complex_simplex_range complex_simplex_range() const { + Complex_simplex_iterator end(this, true); + if (empty()) { + return Complex_simplex_range(end, end); + } else { + Complex_simplex_iterator begin(this); + return Complex_simplex_range(begin, end); + } + } + + //@} + + /** @name Blockers iterators + */ + //@{ + private: + /** + * @brief Iterator over the blockers adjacent to a vertex + */ + typedef Blocker_iterator_around_vertex_internal< + typename std::multimap<Vertex_handle, Simplex *>::iterator, + Blocker_handle> + Complex_blocker_around_vertex_iterator; + + /** + * @brief Iterator over (constant) blockers adjacent to a vertex + */ + typedef Blocker_iterator_around_vertex_internal< + typename std::multimap<Vertex_handle, Simplex *>::const_iterator, + const Blocker_handle> + Const_complex_blocker_around_vertex_iterator; + + typedef boost::iterator_range <Complex_blocker_around_vertex_iterator> Complex_blocker_around_vertex_range; + typedef boost::iterator_range <Const_complex_blocker_around_vertex_iterator> + Const_complex_blocker_around_vertex_range; + + public: + /** + * @brief Returns a range of the blockers of the complex passing through a vertex + */ + Complex_blocker_around_vertex_range blocker_range(Vertex_handle v) { + auto begin = Complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); + auto end = Complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); + return Complex_blocker_around_vertex_range(begin, end); + } + + /** + * @brief Returns a range of the blockers of the complex passing through a vertex + */ + Const_complex_blocker_around_vertex_range const_blocker_range(Vertex_handle v) const { + auto begin = Const_complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v)); + auto end = Const_complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v)); + return Const_complex_blocker_around_vertex_range(begin, end); + } + + private: + /** + * @brief Iterator over the blockers. + */ + typedef Blocker_iterator_internal< + typename std::multimap<Vertex_handle, Simplex *>::iterator, + Blocker_handle> + Complex_blocker_iterator; + + /** + * @brief Iterator over the (constant) blockers. + */ + typedef Blocker_iterator_internal< + typename std::multimap<Vertex_handle, Simplex *>::const_iterator, + const Blocker_handle> + Const_complex_blocker_iterator; + + typedef boost::iterator_range <Complex_blocker_iterator> Complex_blocker_range; + typedef boost::iterator_range <Const_complex_blocker_iterator> Const_complex_blocker_range; + + public: + /** + * @brief Returns a range of the blockers of the complex + */ + Complex_blocker_range blocker_range() { + auto begin = Complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end()); + auto end = Complex_blocker_iterator(blocker_map_.end(), blocker_map_.end()); + return Complex_blocker_range(begin, end); + } + + /** + * @brief Returns a range of the blockers of the complex + */ + Const_complex_blocker_range const_blocker_range() const { + auto begin = Const_complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end()); + auto end = Const_complex_blocker_iterator(blocker_map_.end(), blocker_map_.end()); + return Const_complex_blocker_range(begin, end); + } + + //@} + + ///////////////////////////////////////////////////////////////////////////// + /** @name Print and IO methods + */ + //@{ + public: + std::string to_string() const { + std::ostringstream stream; + stream << num_vertices() << " vertices:\n" << vertices_to_string() << std::endl; + stream << num_edges() << " edges:\n" << edges_to_string() << std::endl; + stream << num_blockers() << " blockers:\n" << blockers_to_string() << std::endl; + return stream.str(); + } + + std::string vertices_to_string() const { + std::ostringstream stream; + for (auto vertex : vertex_range()) { + stream << "{" << (*this)[vertex].get_id() << "} "; + } + stream << std::endl; + return stream.str(); + } + + std::string edges_to_string() const { + std::ostringstream stream; + for (auto edge : edge_range()) + stream << "{" << (*this)[edge].first() << "," << (*this)[edge].second() << "} "; + stream << std::endl; + return stream.str(); + } + + std::string blockers_to_string() const { + std::ostringstream stream; + + for (auto b : const_blocker_range()) + stream << *b << std::endl; + return stream.str(); + } + //@} +}; + +/** + * build a simplicial complex from a collection + * of top faces. + * return the total number of simplices + */ +template<typename Complex, typename SimplexHandleIterator> +Complex make_complex_from_top_faces(SimplexHandleIterator simplices_begin, SimplexHandleIterator simplices_end, + bool is_flag_complex = false) { + // TODO(DS): use add_simplex instead! should be more efficient and more elegant :) + typedef typename Complex::Simplex Simplex; + std::vector<Simplex> simplices; + for (auto top_face = simplices_begin; top_face != simplices_end; ++top_face) { + auto subfaces_topface = subfaces(*top_face); + simplices.insert(simplices.end(), subfaces_topface.begin(), subfaces_topface.end()); + } + return Complex(simplices.begin(), simplices.end(), is_flag_complex); +} + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#include "Skeleton_blocker_simplifiable_complex.h" + +#endif // SKELETON_BLOCKER_COMPLEX_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h new file mode 100644 index 00000000..b8f75e0f --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_geometric_complex.h @@ -0,0 +1,215 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_ +#define SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_ + +#include <gudhi/Skeleton_blocker_complex.h> +#include <gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h> +#include <gudhi/Debug_utils.h> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + * @brief Class that represents a geometric complex that can be simplified. + * The class allows access to points of vertices. + * @ingroup skbl + */ +template<typename SkeletonBlockerGeometricDS> +class Skeleton_blocker_geometric_complex : +public Skeleton_blocker_complex<SkeletonBlockerGeometricDS> { + public: + typedef typename SkeletonBlockerGeometricDS::GT GT; + + typedef Skeleton_blocker_complex<SkeletonBlockerGeometricDS> SimplifiableSkeletonblocker; + + typedef typename SimplifiableSkeletonblocker::Vertex_handle Vertex_handle; + typedef typename SimplifiableSkeletonblocker::Root_vertex_handle Root_vertex_handle; + typedef typename SimplifiableSkeletonblocker::Edge_handle Edge_handle; + typedef typename SimplifiableSkeletonblocker::Simplex Simplex; + + typedef typename SimplifiableSkeletonblocker::Graph_vertex Graph_vertex; + + typedef typename SkeletonBlockerGeometricDS::Point Point; + + Skeleton_blocker_geometric_complex() { } + + /** + * constructor given a list of points + */ + template<typename PointIterator> + explicit Skeleton_blocker_geometric_complex(int num_vertices, PointIterator begin, PointIterator end) { + for (auto point = begin; point != end; ++point) + add_vertex(*point); + } + + /** + * @brief Constructor with a list of simplices. + * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be + * computed or not). + */ + template<typename SimpleHandleOutputIterator, typename PointIterator> + Skeleton_blocker_geometric_complex( + SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, + PointIterator points_begin, PointIterator points_end, + bool is_flag_complex = false) + : Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin, simplex_end, is_flag_complex) { + unsigned current = 0; + for (auto point = points_begin; point != points_end; ++point) + (*this)[Vertex_handle(current++)].point() = Point(point->begin(), point->end()); + } + + /** + * @brief Constructor with a list of simplices. + * Points of every vertex are the point constructed with default constructor. + * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not). + */ + template<typename SimpleHandleOutputIterator> + Skeleton_blocker_geometric_complex( + SimpleHandleOutputIterator simplex_begin, SimpleHandleOutputIterator simplex_end, + bool is_flag_complex = false) + : Skeleton_blocker_complex<SkeletonBlockerGeometricDS>(simplex_begin, simplex_end, is_flag_complex) { } + + /** + * @brief Add a vertex to the complex with a default constructed associated point. + */ + Vertex_handle add_vertex() { + return SimplifiableSkeletonblocker::add_vertex(); + } + + /** + * @brief Add a vertex to the complex with its associated point. + */ + Vertex_handle add_vertex(const Point& point) { + Vertex_handle ad = SimplifiableSkeletonblocker::add_vertex(); + (*this)[ad].point() = point; + return ad; + } + + /** + * @brief Returns the Point associated to the vertex v. + */ + const Point& point(Vertex_handle v) const { + assert(this->contains_vertex(v)); + return (*this)[v].point(); + } + + /** + * @brief Returns the Point associated to the vertex v. + */ + Point& point(Vertex_handle v) { + assert(this->contains_vertex(v)); + return (*this)[v].point(); + } + + const Point& point(Root_vertex_handle global_v) const { + Vertex_handle local_v((*this)[global_v]); + assert(this->contains_vertex(local_v)); + return (*this)[local_v].point(); + } + + Point& point(Root_vertex_handle global_v) { + Vertex_handle local_v((*this)[global_v]); + assert(this->contains_vertex(local_v)); + return (*this)[local_v].point(); + } + + typedef Skeleton_blocker_link_complex<Skeleton_blocker_geometric_complex> Geometric_link; + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Geometric_link link(Vertex_handle v) const { + Geometric_link link(*this, Simplex(v)); + // we now add the point info + add_points_to_link(link); + return link; + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Geometric_link link(const Simplex& simplex) const { + Geometric_link link(*this, simplex); + // we now add the point info + add_points_to_link(link); + return link; + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Geometric_link link(Edge_handle edge) const { + Geometric_link link(*this, edge); + // we now add the point info + add_points_to_link(link); + return link; + } + + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerGeometricDS>> Abstract_link; + + /** + * Constructs the abstract link of v (without points coordinates). + */ + Abstract_link abstract_link(Vertex_handle v) const { + return Abstract_link(*this, Simplex(v)); + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Abstract_link abstract_link(const Simplex& simplex) const { + return Abstract_link(*this, simplex); + } + + /** + * Constructs the link of 'simplex' with points coordinates. + */ + Abstract_link abstract_link(Edge_handle edge) const { + return Abstract_link(*this, edge); + } + + private: + void add_points_to_link(Geometric_link& link) const { + for (Vertex_handle v : link.vertex_range()) { + Root_vertex_handle v_root(link.get_id(v)); + link.point(v) = (*this).point(v_root); + } + } +}; + + +template<typename SkeletonBlockerGeometricComplex, typename SimpleHandleOutputIterator, typename PointIterator> +SkeletonBlockerGeometricComplex make_complex_from_top_faces( + SimpleHandleOutputIterator simplex_begin, + SimpleHandleOutputIterator simplex_end, + PointIterator points_begin, + PointIterator points_end, + bool is_flag_complex = false) { + typedef SkeletonBlockerGeometricComplex SBGC; + SkeletonBlockerGeometricComplex complex; + unsigned current = 0; + complex = + make_complex_from_top_faces<SBGC>(simplex_begin, simplex_end, is_flag_complex); + for (auto point = points_begin; point != points_end; ++point) + // complex.point(Vertex_handle(current++)) = Point(point->begin(),point->end()); + complex.point(typename SBGC::Vertex_handle(current++)) = typename SBGC::Point(*point); + return complex; +} + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_GEOMETRIC_COMPLEX_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_link_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_link_complex.h new file mode 100644 index 00000000..a2637da3 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_link_complex.h @@ -0,0 +1,287 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_LINK_COMPLEX_H_ +#define SKELETON_BLOCKER_LINK_COMPLEX_H_ + +#include <gudhi/Skeleton_blocker_complex.h> +#include <gudhi/Debug_utils.h> + +namespace Gudhi { + +namespace skeleton_blocker { + +template<class ComplexType> class Skeleton_blocker_sub_complex; + +/** + * \brief Class representing the link of a simplicial complex encoded by a skeleton/blockers pair. + * It inherits from Skeleton_blocker_sub_complex because such complex is a sub complex of a + * root complex. + * \ingroup skbl + */ +template<typename ComplexType> +class Skeleton_blocker_link_complex : public Skeleton_blocker_sub_complex< +ComplexType> { + template<typename T> friend class Skeleton_blocker_link_superior; + typedef typename ComplexType::Edge_handle Edge_handle; + + typedef typename ComplexType::boost_vertex_handle boost_vertex_handle; + + private: + bool only_superior_vertices_; + + public: + typedef typename ComplexType::Vertex_handle Vertex_handle; + typedef typename ComplexType::Root_vertex_handle Root_vertex_handle; + + typedef typename ComplexType::Simplex Simplex; + typedef typename ComplexType::Root_simplex_handle Root_simplex_handle; + + typedef typename ComplexType::Blocker_handle Blocker_handle; + + typedef typename ComplexType::Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_handle_iterator; + + explicit Skeleton_blocker_link_complex(bool only_superior_vertices = false) + : only_superior_vertices_(only_superior_vertices) { } + + /** + * If the parameter only_superior_vertices is true, + * only vertices greater than the one of alpha are added. + * Only vertices are computed if only_vertices is true. + */ + Skeleton_blocker_link_complex(const ComplexType & parent_complex, + const Simplex& alpha_parent_adress, + bool only_superior_vertices = false, + bool only_vertices = false) + : only_superior_vertices_(only_superior_vertices) { + if (!alpha_parent_adress.empty()) + build_link(parent_complex, alpha_parent_adress, only_vertices); + } + + /** + * If the parameter only_superior_vertices is true, + * only vertices greater than the one of the vertex are added. + */ + Skeleton_blocker_link_complex(const ComplexType & parent_complex, + Vertex_handle a_parent_adress, + bool only_superior_vertices = false) + : only_superior_vertices_(only_superior_vertices) { + Simplex alpha_simplex(a_parent_adress); + build_link(parent_complex, alpha_simplex); + } + + /** + * If the parameter only_superior_vertices is true, + * only vertices greater than the one of the edge are added. + */ + Skeleton_blocker_link_complex(const ComplexType & parent_complex, + Edge_handle edge, bool only_superior_vertices = + false) + : only_superior_vertices_(only_superior_vertices) { + Simplex alpha_simplex(parent_complex.first_vertex(edge), + parent_complex.second_vertex(edge)); + build_link(parent_complex, alpha_simplex); + } + + protected: + /** + * @brief compute vertices of the link. + * If the boolean only_superior_vertices is true, then only the vertices + * are greater than vertices of alpha_parent_adress are added. + */ + void compute_link_vertices(const ComplexType & parent_complex, + const Simplex& alpha_parent_adress, + bool only_superior_vertices, + bool is_alpha_blocker = false) { + if (alpha_parent_adress.dimension() == 0) { + // for a vertex we know exactly the number of vertices of the link (and the size of the corresponding vector) + // thus we call a specific function that will reserve a vector with appropriate size + this->compute_link_vertices(parent_complex, + alpha_parent_adress.first_vertex(), + only_superior_vertices_); + } else { + // we compute the intersection of neighbors of alpha and store it in link_vertices + Simplex link_vertices_parent; + parent_complex.add_neighbours(alpha_parent_adress, link_vertices_parent, + only_superior_vertices); + // For all vertex 'v' in this intersection, we go through all its adjacent blockers. + // If one blocker minus 'v' is included in alpha then the vertex is not in the link complex. + for (auto v_parent : link_vertices_parent) { + bool new_vertex = true; + for (auto beta : parent_complex.const_blocker_range(v_parent)) { + if (!is_alpha_blocker || *beta != alpha_parent_adress) { + new_vertex = !(alpha_parent_adress.contains_difference(*beta, + v_parent)); + if (!new_vertex) + break; + } + } + if (new_vertex) + this->add_vertex(parent_complex.get_id(v_parent)); + } + } + } + + /** + * @brief compute vertices of the link. + * If the boolean only_superior_vertices is true, then only the vertices + * are greater than vertices of alpha_parent_adress are added. + */ + void compute_link_vertices(const ComplexType & parent_complex, + Vertex_handle alpha_parent_adress, + bool only_superior_vertices) { + // for a vertex we know exactly the number of vertices of the link (and the size of the corresponding vector + this->skeleton.m_vertices.reserve( + parent_complex.degree(alpha_parent_adress)); + + // For all vertex 'v' in this intersection, we go through all its adjacent blockers. + // If one blocker minus 'v' is included in alpha then the vertex is not in the link complex. + for (auto v_parent : parent_complex.vertex_range(alpha_parent_adress)) { + if (!only_superior_vertices + || v_parent.vertex > alpha_parent_adress.vertex) + this->add_vertex(parent_complex.get_id(v_parent)); + } + } + + void compute_link_edges(const ComplexType & parent_complex, + const Simplex& alpha_parent_adress, + bool is_alpha_blocker = false) { + if (this->num_vertices() <= 1) + return; + + for (auto x_link = this->vertex_range().begin(); + x_link != this->vertex_range().end(); ++x_link) { + for (auto y_link = x_link; ++y_link != this->vertex_range().end();) { + Vertex_handle x_parent = *parent_complex.get_address( + this->get_id(*x_link)); + Vertex_handle y_parent = *parent_complex.get_address( + this->get_id(*y_link)); + if (parent_complex.contains_edge(x_parent, y_parent)) { + // we check that there is no blocker subset of alpha passing trough x and y + bool new_edge = true; + for (auto blocker_parent : parent_complex.const_blocker_range( + x_parent)) { + if (!is_alpha_blocker || *blocker_parent != alpha_parent_adress) { + if (blocker_parent->contains(y_parent)) { + new_edge = !(alpha_parent_adress.contains_difference( + *blocker_parent, x_parent, y_parent)); + if (!new_edge) + break; + } + } + } + if (new_edge) + this->add_edge_without_blockers(*x_link, *y_link); + } + } + } + } + + /** + * @brief : Given an address in the current complex, it returns the + * corresponding address in 'other_complex'. + * It assumes that other_complex have a vertex 'this.get_id(address)' + */ + boost::optional<Vertex_handle> give_equivalent_vertex(const ComplexType & other_complex, + Vertex_handle address) const { + Root_vertex_handle id((*this)[address].get_id()); + return other_complex.get_address(id); + } + + /* + * compute the blockers of the link if is_alpha_blocker is false. + * Otherwise, alpha is a blocker, and the link is computed in the complex where + * the blocker is anticollapsed. + */ + void compute_link_blockers(const ComplexType & parent_complex, + const Simplex& alpha_parent, + bool is_alpha_blocker = false) { + for (auto x_link : this->vertex_range()) { + Vertex_handle x_parent = *this->give_equivalent_vertex(parent_complex, + x_link); + + for (auto blocker_parent : parent_complex.const_blocker_range(x_parent)) { + if (!is_alpha_blocker || *blocker_parent != alpha_parent) { + Simplex sigma_parent(*blocker_parent); + + sigma_parent.difference(alpha_parent); + + if (sigma_parent.dimension() >= 2 + && sigma_parent.first_vertex() == x_parent) { + Root_simplex_handle sigma_id(parent_complex.get_id(sigma_parent)); + auto sigma_link = this->get_simplex_address(sigma_id); + // ie if the vertices of sigma are vertices of the link + if (sigma_link) { + bool is_new_blocker = true; + for (auto a : alpha_parent) { + for (auto eta_parent : parent_complex.const_blocker_range(a)) { + if (!is_alpha_blocker || *eta_parent != alpha_parent) { + Simplex eta_minus_alpha(*eta_parent); + eta_minus_alpha.difference(alpha_parent); + if (eta_minus_alpha != sigma_parent + && sigma_parent.contains_difference(*eta_parent, + alpha_parent)) { + is_new_blocker = false; + break; + } + } + } + if (!is_new_blocker) + break; + } + if (is_new_blocker) + this->add_blocker(new Simplex(*sigma_link)); + } + } + } + } + } + } + + public: + /** + * @brief compute vertices, edges and blockers of the link. + * @details If the boolean only_superior_vertices is true, then the link is computed only + * with vertices that are greater than vertices of alpha_parent_adress. + */ + void build_link(const ComplexType & parent_complex, + const Simplex& alpha_parent_adress, + bool is_alpha_blocker = false, + bool only_vertices = false) { + assert(is_alpha_blocker || parent_complex.contains(alpha_parent_adress)); + compute_link_vertices(parent_complex, alpha_parent_adress, only_superior_vertices_); + if (!only_vertices) { + compute_link_edges(parent_complex, alpha_parent_adress, is_alpha_blocker); + compute_link_blockers(parent_complex, alpha_parent_adress, is_alpha_blocker); + } + } + + /** + * @brief build the link of a blocker which is the link + * of the blocker's simplex if this simplex had been + * removed from the blockers of the complex. + */ + friend void build_link_of_blocker(const ComplexType & parent_complex, + Simplex& blocker, + Skeleton_blocker_link_complex & result) { + assert(blocker.dimension() >= 2); + assert(parent_complex.contains_blocker(blocker)); + result.clear(); + result.build_link(parent_complex, blocker, true); + } +}; + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_LINK_COMPLEX_H_ diff --git a/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h new file mode 100644 index 00000000..404f04f9 --- /dev/null +++ b/src/Skeleton_blocker/include/gudhi/Skeleton_blocker_simplifiable_complex.h @@ -0,0 +1,455 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_SIMPLIFIABLE_COMPLEX_H_ +#define SKELETON_BLOCKER_SIMPLIFIABLE_COMPLEX_H_ + +#include <gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h> + +#include <list> +#include <vector> +#include <set> + +namespace Gudhi { + +namespace skeleton_blocker { + +/** + * Returns true if the blocker 'sigma' is popable. + * To define popable, let us call 'L' the complex that + * consists in the current complex without the blocker 'sigma'. + * A blocker 'sigma' is then "popable" if the link of 'sigma' + * in L is reducible. + * + */ +template<typename SkeletonBlockerDS> +bool Skeleton_blocker_complex<SkeletonBlockerDS>::is_popable_blocker(Blocker_handle sigma) const { + assert(this->contains_blocker(*sigma)); + Skeleton_blocker_link_complex<Skeleton_blocker_complex> link_blocker_sigma; + build_link_of_blocker(*this, *sigma, link_blocker_sigma); + return link_blocker_sigma.is_contractible() == CONTRACTIBLE; +} + +/** + * Removes all the popable blockers of the complex and delete them. + * @returns the number of popable blockers deleted + */ +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers() { + std::list<Vertex_handle> vertex_to_check; + for (auto v : this->vertex_range()) + vertex_to_check.push_front(v); + + while (!vertex_to_check.empty()) { + Vertex_handle v = vertex_to_check.front(); + vertex_to_check.pop_front(); + + bool blocker_popable_found = true; + while (blocker_popable_found) { + blocker_popable_found = false; + for (auto block : this->blocker_range(v)) { + if (this->is_popable_blocker(block)) { + for (Vertex_handle nv : *block) + if (nv != v) vertex_to_check.push_back(nv); + this->delete_blocker(block); + blocker_popable_found = true; + break; + } + } + } + } +} + +/** + * Removes all the popable blockers of the complex passing through v and delete them. + */ +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_popable_blockers(Vertex_handle v) { + bool blocker_popable_found = true; + while (blocker_popable_found) { + blocker_popable_found = false; + for (auto block : this->blocker_range(v)) { + if (is_popable_blocker(block)) { + this->delete_blocker(block); + blocker_popable_found = true; + } + } + } +} + +/** + * @brief Removes all the popable blockers of the complex passing through v and delete them. + * Also remove popable blockers in the neighborhood if they became popable. + * + */ +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_all_popable_blockers(Vertex_handle v) { + std::list<Vertex_handle> vertex_to_check; + vertex_to_check.push_front(v); + + while (!vertex_to_check.empty()) { + Vertex_handle v = vertex_to_check.front(); + vertex_to_check.pop_front(); + + bool blocker_popable_found = true; + while (blocker_popable_found) { + blocker_popable_found = false; + for (auto block : this->blocker_range(v)) { + if (this->is_popable_blocker(block)) { + for (Vertex_handle nv : *block) + if (nv != v) vertex_to_check.push_back(nv); + this->delete_blocker(block); + blocker_popable_found = true; + break; + } + } + } + } +} + +/** + * Remove the star of the vertice 'v' + */ +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle v) { + // we remove the blockers that are not consistent anymore + update_blockers_after_remove_star_of_vertex_or_edge(Simplex(v)); + while (this->degree(v) > 0) { + Vertex_handle w(* (adjacent_vertices(v.vertex, this->skeleton).first)); + this->remove_edge(v, w); + } + this->remove_vertex(v); +} + +/** + * after removing the star of a simplex, blockers sigma that contains this simplex must be removed. + * Furthermore, all simplices tau of the form sigma \setminus simplex_to_be_removed must be added + * whenever the dimension of tau is at least 2. + */ +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::update_blockers_after_remove_star_of_vertex_or_edge(const Simplex& simplex_to_be_removed) { + std::list <Blocker_handle> blockers_to_update; + if (simplex_to_be_removed.empty()) return; + + auto v0 = simplex_to_be_removed.first_vertex(); + for (auto blocker : this->blocker_range(v0)) { + if (blocker->contains(simplex_to_be_removed)) + blockers_to_update.push_back(blocker); + } + + for (auto blocker_to_update : blockers_to_update) { + Simplex sub_blocker_to_be_added; + bool sub_blocker_need_to_be_added = + (blocker_to_update->dimension() - simplex_to_be_removed.dimension()) >= 2; + if (sub_blocker_need_to_be_added) { + sub_blocker_to_be_added = *blocker_to_update; + sub_blocker_to_be_added.difference(simplex_to_be_removed); + } + this->delete_blocker(blocker_to_update); + if (sub_blocker_need_to_be_added) + this->add_blocker(sub_blocker_to_be_added); + } +} + +/** + * Remove the star of the edge connecting vertices a and b. + * @returns the number of blocker that have been removed + */ +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Vertex_handle a, Vertex_handle b) { + update_blockers_after_remove_star_of_vertex_or_edge(Simplex(a, b)); + // we remove the edge + this->remove_edge(a, b); +} + +/** + * Remove the star of the edge 'e'. + */ +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(Edge_handle e) { + return remove_star(this->first_vertex(e), this->second_vertex(e)); +} + +/** + * Remove the star of the simplex 'sigma' which needs to belong to the complex + */ +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_star(const Simplex& sigma) { + assert(this->contains(sigma)); + if (sigma.dimension() == 0) { + remove_star(sigma.first_vertex()); + } else if (sigma.dimension() == 1) { + remove_star(sigma.first_vertex(), sigma.last_vertex()); + } else { + remove_blocker_containing_simplex(sigma); + this->add_blocker(sigma); + } +} + +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::add_simplex(const Simplex& sigma) { + // to add a simplex s, all blockers included in s are first removed + // and then all simplex in the coboundary of s are added as blockers + assert(!this->contains(sigma)); + assert(sigma.dimension() > 1); + if (!contains_vertices(sigma)) { + std::cerr << "add_simplex: Some vertices were not present in the complex, adding them" << std::endl; + size_t num_vertices_to_add = sigma.last_vertex() - this->num_vertices() + 1; + for (size_t i = 0; i < num_vertices_to_add; ++i) + this->add_vertex(); + } + assert(contains_vertices(sigma)); + if (!contains_edges(sigma)) + add_edge(sigma); + remove_blocker_include_in_simplex(sigma); + add_blockers_after_simplex_insertion(sigma); +} + +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::add_blockers_after_simplex_insertion(Simplex sigma) { + if (sigma.dimension() < 1) return; + + for (auto s : coboundary_range(sigma)) { + this->add_blocker(s); + } +} + +/** + * remove all blockers that contains sigma + */ +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_containing_simplex(const Simplex& sigma) { + std::vector <Blocker_handle> blockers_to_remove; + for (auto blocker : this->blocker_range(sigma.first_vertex())) { + if (blocker->contains(sigma)) + blockers_to_remove.push_back(blocker); + } + for (auto blocker_to_update : blockers_to_remove) + this->delete_blocker(blocker_to_update); +} + +/** + * remove all blockers that contains sigma + */ +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::remove_blocker_include_in_simplex(const Simplex& sigma) { + // TODO(DS): write efficiently by using only superior blockers + // eg for all s, check blockers whose vertices are all greater than s + std::set <Blocker_handle> blockers_to_remove; + for (auto s : sigma) { + for (auto blocker : this->blocker_range(s)) { + if (sigma.contains(*blocker)) + blockers_to_remove.insert(blocker); + } + } + for (auto blocker_to_update : blockers_to_remove) { + auto s = *blocker_to_update; + this->delete_blocker(blocker_to_update); + // now if there is a vertex v in the link of s + // and v is not included in sigma then v.s is a blocker + // (all faces of v.s are there since v belongs to the link of s) + for (const auto& b : coboundary_range(s)) + if (!sigma.contains(b)) + this->add_blocker(b); + } +} + +/** + * Compute simplices beta such that a.beta is an order 0 blocker + * that may be used to construct a new blocker after contracting ab. + * It requires that the link condition is satisfied. + */ +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::tip_blockers(Vertex_handle a, Vertex_handle b, + std::vector<Simplex> & buffer) const { + for (auto const & blocker : this->const_blocker_range(a)) { + Simplex beta = (*blocker); + beta.remove_vertex(a); + buffer.push_back(beta); + } + + Simplex n; + this->add_neighbours(b, n); + this->remove_neighbours(a, n); + n.remove_vertex(a); + + + for (Vertex_handle y : n) { + Simplex beta; + beta.add_vertex(y); + buffer.push_back(beta); + } +} + +/** + * @brief "Replace" the edge 'bx' by the edge 'ax'. + * Assume that the edge 'bx' was present whereas 'ax' was not. + * Precisely, it does not replace edges, but remove 'bx' and then add 'ax'. + * The visitor 'on_swaped_edge' is called just after edge 'ax' had been added + * and just before edge 'bx' had been removed. That way, it can + * eventually access to information of 'ax'. + */ +template<typename SkeletonBlockerDS> +void +Skeleton_blocker_complex<SkeletonBlockerDS>::swap_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x) { + this->add_edge_without_blockers(a, x); + if (this->visitor) this->visitor->on_swaped_edge(a, b, x); + this->remove_edge(b, x); +} + +template<typename SkeletonBlockerDS> +void +Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_vertex(Vertex_handle v) { + std::list <Blocker_handle> blockers_to_delete; + for (auto blocker : this->blocker_range(v)) { + blockers_to_delete.push_back(blocker); + } + while (!blockers_to_delete.empty()) { + this->remove_blocker(blockers_to_delete.back()); + blockers_to_delete.pop_back(); + } +} + +/** + * @brief removes all blockers passing through the edge 'ab' + */ +template<typename SkeletonBlockerDS> +void +Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_edge(Vertex_handle a, Vertex_handle b) { + std::list<Blocker_handle> blocker_to_delete; + for (auto blocker : this->blocker_range(a)) + if (blocker->contains(b)) blocker_to_delete.push_back(blocker); + while (!blocker_to_delete.empty()) { + this->delete_blocker(blocker_to_delete.back()); + blocker_to_delete.pop_back(); + } +} + +/** + * Contracts the edge connecting vertices a and b. + * @remark If the link condition Link(ab) = Link(a) inter Link(b) is not satisfied, + * it removes first all blockers passing through 'ab' + */ +template<typename SkeletonBlockerDS> +void +Skeleton_blocker_complex<SkeletonBlockerDS>::contract_edge(Vertex_handle a, Vertex_handle b) { + assert(this->contains_vertex(a)); + assert(this->contains_vertex(b)); + + if (this->contains_edge(a, b)) + this->add_edge_without_blockers(a, b); + + // if some blockers passes through 'ab', we need to remove them. + if (!link_condition(a, b)) + delete_blockers_around_edge(a, b); + + std::set<Simplex> blockers_to_add; + + get_blockers_to_be_added_after_contraction(a, b, blockers_to_add); + + delete_blockers_around_vertices(a, b); + + update_edges_after_contraction(a, b); + + this->remove_vertex(b); + + notify_changed_edges(a); + + for (auto block : blockers_to_add) + this->add_blocker(block); + + assert(this->contains_vertex(a)); + assert(!this->contains_vertex(b)); +} + +template<typename SkeletonBlockerDS> +void Skeleton_blocker_complex<SkeletonBlockerDS>::get_blockers_to_be_added_after_contraction(Vertex_handle a, + Vertex_handle b, + std::set<Simplex>& blockers_to_add) { + blockers_to_add.clear(); + + typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > LinkComplexType; + + LinkComplexType link_a(*this, a); + LinkComplexType link_b(*this, b); + + std::vector<Simplex> vector_alpha, vector_beta; + + tip_blockers(a, b, vector_alpha); + tip_blockers(b, a, vector_beta); + + for (auto alpha = vector_alpha.begin(); alpha != vector_alpha.end(); ++alpha) { + for (auto beta = vector_beta.begin(); beta != vector_beta.end(); ++beta) { + Simplex sigma = *alpha; + sigma.union_vertices(*beta); + Root_simplex_handle sigma_id = this->get_id(sigma); + if (this->contains(sigma) && + proper_faces_in_union<Skeleton_blocker_complex < SkeletonBlockerDS >> (sigma_id, link_a, link_b)) { + // Blocker_handle blocker = new Simplex(sigma); + sigma.add_vertex(a); + blockers_to_add.insert(sigma); + } + } + } +} + +/** + * delete all blockers that passes through a or b + */ +template<typename SkeletonBlockerDS> +void +Skeleton_blocker_complex<SkeletonBlockerDS>::delete_blockers_around_vertices(Vertex_handle a, Vertex_handle b) { + std::vector<Blocker_handle> blocker_to_delete; + for (auto bl : this->blocker_range(a)) + blocker_to_delete.push_back(bl); + for (auto bl : this->blocker_range(b)) + blocker_to_delete.push_back(bl); + while (!blocker_to_delete.empty()) { + this->delete_blocker(blocker_to_delete.back()); + blocker_to_delete.pop_back(); + } +} + +template<typename SkeletonBlockerDS> +void +Skeleton_blocker_complex<SkeletonBlockerDS>::update_edges_after_contraction(Vertex_handle a, Vertex_handle b) { + // We update the set of edges + this->remove_edge(a, b); + + // For all edges {b,x} incident to b, + // we remove {b,x} and add {a,x} if not already there. + while (this->degree(b) > 0) { + Vertex_handle x(*(adjacent_vertices(b.vertex, this->skeleton).first)); + if (!this->contains_edge(a, x)) + // we 'replace' the edge 'bx' by the edge 'ax' + this->swap_edge(a, b, x); + else + this->remove_edge(b, x); + } +} + +template<typename SkeletonBlockerDS> +void +Skeleton_blocker_complex<SkeletonBlockerDS>::notify_changed_edges(Vertex_handle a) { + // We notify the visitor that all edges incident to 'a' had changed + boost_adjacency_iterator v, v_end; + + for (tie(v, v_end) = adjacent_vertices(a.vertex, this->skeleton); v != v_end; ++v) + if (this->visitor) this->visitor->on_changed_edge(a, Vertex_handle(*v)); +} + + +} // namespace skeleton_blocker + +namespace skbl = skeleton_blocker; + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_SIMPLIFIABLE_COMPLEX_H_ |